A coarse-grained NADH redox model enables inference of subcellular metabolic fluxes from fluorescence lifetime imaging

Mitochondrial metabolism is of central importance to diverse aspects of cell and developmental biology. Defects in mitochondria are associated with many diseases, including cancer, neuropathology, and infertility. Our understanding of mitochondrial metabolism in situ and dysfunction in diseases are limited by the lack of techniques to measure mitochondrial metabolic fluxes with sufficient spatiotemporal resolution. Herein, we developed a new method to infer mitochondrial metabolic fluxes in living cells with subcellular resolution from fluorescence lifetime imaging of NADH. This result is based on the use of a generic coarse-grained NADH redox model. We tested the model in mouse oocytes and human tissue culture cells subject to a wide variety of perturbations by comparing predicted fluxes through the electron transport chain (ETC) to direct measurements of oxygen consumption rate. Interpreting the fluorescence lifetime imaging microscopy measurements of NADH using this model, we discovered a homeostasis of ETC flux in mouse oocytes: perturbations of nutrient supply and energy demand of the cell do not change ETC flux despite significantly impacting NADH metabolic state. Furthermore, we observed a subcellular spatial gradient of ETC flux in mouse oocytes and found that this gradient is primarily a result of a spatially heterogeneous mitochondrial proton leak. We concluded from these observations that ETC flux in mouse oocytes is not controlled by energy demand or supply, but by the intrinsic rates of mitochondrial respiration.


Introduction
Cells transduce energy from the environment to power cellular processes. Decades of extensive research have produced a remarkable body of detailed information about the biochemistry of mitochondrial energy metabolism (Salway, 2017). In brief, metabolites, such as pyruvate, are transported into mitochondria, where they are broken down and their products enter the tricarboxylic acid cycle (TCA). The TCA is composed of a number of chemical reactions, which ultimately reduces NAD + to NADH. NADH and oxygen are then utilized by the electron transport chain (ETC) to pump hydrogen ions across the mitochondrial membrane. ATP synthase uses this proton gradient to power the synthesis of ATP from ADP (Mitchell, 1961). The activities of mitochondrial energy metabolism are characterized by the fluxes through these pathways: that is, the number of molecules turned over per unit time (Stephanopoulos, 1999). However, despite the wealth of knowledge concerning mitochondrial biochemistry, the spatiotemporal dynamics of cellular energy usage remains elusive and it is still unclear how cells partition energy across different cellular processes (Dumollard et al., 2007;Van Blerkom, 2011;Yellen, 2018;Yang et al., 2021) and how energy metabolism is misregulated in diseases (Brand and Nicholls, 2011;Lin and Beal, 2006;Wallace, 2012;Bratic and Larsson, 2013;Lowell and Shulman, 2005;Mick et al., 2020). Metabolic heterogeneities, between and within individual cells, are believed to be widespread, but remain poorly characterized (Takhaveev and Heinemann, 2018;Aryaman et al., 2018). Mitochondria have been observed to associate with the cytoskeleton (Lawrence et al., 2016), spindle (Wang et al., 2020), and endoplasmic reticulum (Dumollard et al., 2004) and display subcellular heterogeneities in mtDNA sequence (Morris et al., 2017) and mitochondrial membrane potential (Smiley et al., 1991). These observations suggest the potential existence of subcellular patterning of mitochondrial metabolic fluxes that could be critical in processes such as oocyte maturation (Yu et al., 2010) and embryo development (Sanchez et al., 2019). The limitations of current techniques for measuring mitochondrial metabolic fluxes with sufficient spatiotemporal resolution present a major challenge. In particular, there is a lack of techniques to measure mitochondrial metabolic fluxes with single cell and subcellular resolution.
Bulk biochemical techniques for measuring metabolic fluxes, such as oxygen consumption and nutrient uptake rates (Ferrick et al., 2008;Houghton et al., 1996;Lopes et al., 2005), and isotope tracing by mass spectrometry (Wiechert, 2001), require averaging over large populations of cells. Such techniques cannot resolve cellular, or subcellular, metabolic heterogeneity (Takhaveev and Heinemann, 2018;Aryaman et al., 2018). Biochemical approaches for measuring mitochondrial metabolic fluxes, such as mass spectrometry, are also often destructive (Wiechert, 2001;Saks et al., 1998), and thus cannot be used to observe continual changes in fluxes over time. Fluorescence microscopy provides a powerful means to measure cellular and subcellular metabolic heterogeneity continuously and non-destructively, with high spatiotemporal resolution. However, while fluorescent probes can be used to measure mitochondrial membrane potential (Perry et al., 2011) and the concentration of key metabolites (Imamura et al., 2009;Berg et al., 2009;Díaz-García et al., 2017;San Martín et al., 2014), it is not clear how to relate those observables to mitochondrial metabolic fluxes.
NADH is an important cofactor that is involved in many metabolic pathways, including the TCA and ETC in mitochondria. NADH binds with enzymes and acts as an electron carrier that facilitates redox reactions. In the ETC, for example, NADH binds to complex I and donates its electron to ubiquinone and ultimately to oxygen, becoming oxidized to NAD + . Endogenous NADH has long been used to non-invasively probe cellular metabolism because NADH is autofluorescent, while NAD + is not (Heikal, 2010). Fluorescence lifetime imaging microscopy (FLIM) of NADH autofluorescence allows quantitative measurements of the concentration of NADH, the fluorescence lifetimes of NADH, and the fraction of NADH molecules bound to enzymes (Becker, 2012;Becker, 2019;Bird et al., 2005;Skala et al., 2007;Heikal, 2010;Sharick et al., 2018;Sanchez et al., 2018;Sanchez et al., 2019;Ma et al., 2019). It has been observed that the fraction of enzyme-bound NADH and NADH fluorescence lifetimes are correlated with the activity of oxidative phosphorylation, indicating that there is a connection between NADH enzyme-binding and mitochondrial metabolic fluxes (Bird et al., 2005;Skala et al., 2007). The mechanistic basis of this empirical correlation has been unclear.
Here, we developed a generic coarse-grained NADH redox model that enables the inference of ETC flux with subcellular resolution from FLIM measurements. We validated this model in mouse oocytes and human tissue culture cells subject to a wide range of perturbations by comparing predicted ETC fluxes from FLIM to direct measurements of oxygen consumption rate (OCR), and by a self-consistency criterion. Using this method, we discovered that perturbing nutrient supply and energy demand significantly impacts NADH metabolic state but does not change ETC flux. We also discovered a subcellular spatial gradient of ETC flux in mouse oocytes and found that this flux gradient is primarily due to a spatially heterogeneous mitochondrial proton leak. We concluded from these observations that ETC flux in mouse oocytes is not controlled by energy demand or supply, but by the intrinsic rates of mitochondrial respiration. Thus, FLIM of NADH can be used to non-invasively and , and the concentration of enzyme-bound NADH, , in mitochondria as a function of oxygen (lower right). Error bars are standard error of the mean (s.e.m) across individual oocytes. FLIM, fluorescence lifetime imaging microscopy.
The online version of this article includes the following figure supplement(s) for figure 1:   continuously measure mitochondrial ETC fluxes with subcellular resolution and provides novel insights into the spatiotemporal regulation of metabolic fluxes in cells.

Results
Quantifying response of mitochondrial metabolism to changing oxygen levels and metabolic inhibitors using FLIM of NADH We used meiosis II arrested mouse oocytes as a model system. MII oocytes are in a metabolic steadystate, which eases interpretations of metabolic perturbations. ATP synthesis in mouse oocytes occurs primarily through oxidative phosphorylation using pyruvate, without an appreciable contribution from glycolysis (Houghton et al., 1996), providing an excellent system to study mitochondrial metabolism. Mouse oocytes can be cultured in vitro using chemically well-defined media (Biggers and Racowsky, 2002). In our work, we used AKSOM as the culturing media (Summers, 2013). The oocytes can directly take up pyruvate supplied to them or derive it from lactate through the activity of lactate dehydrogenase (LDH) (Lane and Gardner, 2000), and they can remain in a steady-state for hours with constant metabolic fluxes. While NADH and NADPH are difficult to distinguish with fluorescence measurements due to their overlapping fluorescence spectrum, the concentration of NADH in mouse oocytes is 40 times greater than the concentration of NADPH for the whole cell (Bustamante et al., 2017) and potentially even greater in mitochondria (Zhao et al., 2011), so the autofluorescence signal from these cells (particularly from mitochondria) can be safely assumed to result from NADH.
To investigate how FLIM measurements vary with mitochondrial activities, we performed quantitative metabolic perturbations. We first continually varied the concentration of oxygen in the media, from 50±2 μM to 0.26±0.04 μM, over a course of 30 min while imaging NADH autofluorescence of oocytes with FLIM ( Figure 1a, top, black curve; Video 1). NADH is present in both mitochondria and cytoplasm where it is involved in different metabolic pathways. To specifically study the response of NADH in mitochondria, we used a machine learning-based algorithm to segment mitochondria from the NADH intensity images (Berg et al., 2019; Figure 1b and Figure 1-figure supplement 1). We ). (f) Bound NADH concentrations ( ). Error bars represent standard error of the mean (s.e.m) across different oocytes. Student's t-test is performed between parameters before and after the perturbation. *p<0.05, **p<0.01, ***p<0.001. FLIM, fluorescence lifetime imaging microscopy.
The online version of this article includes the following source data for figure 2: Source data 1. Excel spreadsheet of single-oocyte FLIM data used for Figure 2a -f. verified the accuracy of the segmentation with a mitochondrial labeling dye, MitoTracker Red FM, which showed a 80.6±1% (SEM) accuracy of the segmentation (Appendix 1).
Using the segmentation mask, we obtained the intensity of NADH, I , in mitochondria by averaging the photon count over all mitochondrial pixels. The intensity increased with decreasing oxygen concentration (Figure 1a, top, red), as is readily seen from the raw images (Figure 1a, middle). Restoring oxygen to its original level caused a recovery of NADH intensity, indicating that the observed changes are reversible (Figure 1a; Video 1). These observations are consistent with the expectation that NADH concentration will increase at low oxygen levels due to oxygen's role as the electron acceptor in the ETC. In addition to intensity, FLIM can be used to determine the enzyme engagement of NADH by measuring the photon arrival times, from which fluorescence lifetimes can be fitted. Enzyme-bound NADH has a much longer fluorescence lifetime than free NADH (Sharick et al., 2018), allowing bound and free NADH to be separately resolved, but the precise fluorescence lifetimes of NADH depend on a range of factors, including viscosity, pH, and the identity of the enzyme NADH binds to (Sharick et al., 2018;Ghukasyan and Heikal, 2015). To fit NADH fluorescence lifetimes, we grouped all detected photons from mitochondria to form histograms of photon arrival times from NADH autofluorescence for each time point (Figure 1a, lower). We fitted the histograms using a model in which the NADH fluorescence decay, F(τ ) , is described by the sum of two exponentials, where τ l and τs are long and short fluorescence lifetimes, corresponding to enzyme-bound NADH and free NADH, respectively, and f is the fraction of enzyme-bound NADH (Sanchez et al., 2018;Sanchez et al., 2019) (Materials and methods). We repeated the oxygen drop experiments for a total of 68 oocytes. Since the oxygen drop is much slower than the NADH redox reactions (30 min compared to a timescale of seconds), the oxygen perturbation can be safely assumed to be quasistatic, allowing the FLIM measurements to be determined as a function of oxygen levels. We averaged data from all oocytes to obtain a total of four FLIM parameters: mitochondrial NADH intensity, I, long and short fluorescence lifetimes, τ l and τs , and the fraction of enzyme-bound NADH, f . We determined how these parameters varied with oxygen level (Figure 1a and c). All parameters are insensitive to oxygen level until oxygen drops below ~10 μM. This observation is consistent with previous studies that showed mitochondria have a very high apparent affinity for oxygen (Chance and Williams, 1955;Gnaiger et al., 1998).
We next explored the relationship between the measured FLIM parameters and the concentration of NADH. Since bound and free NADH have different fluorescence lifetimes, and hence different molecular brightnesses, the NADH concentration is not generally proportional to NADH intensity. Assuming molecular brightness is proportional to fluorescence lifetime (Lakowicz, 2006), we derived a relation between NADH intensity, fluorescence lifetimes, and concentrations as [ where cs is a calibration factor that relates intensities and concentrations (see Appendix 1). We measured the calibration factor by titrating free NADH in vitro and acquiring FLIM data ( Figure 1 LDH concentration while the sum of free and bound NADH concentration remains a constant and equal to the amount of NADH added to the solution (Figure 1-figure supplement 3). This result demonstrates that Equations 2a,b can be used to measure free and bound NADH concentrations from NADH intensity and lifetimes. It is well established that FLIM can be used to distinguish bound and free NADH in vivo based on the large change of fluorescence lifetime when NADH binds to enzymes (Skala et al., 2007;Heikal, 2010). Even though the exact amount that the lifetime changes depend on the specific enzyme NADH binds to (Sharick et al., 2018), enzyme-bound NADH always has a much longer fluorescence lifetime than free NADH. Therefore, the method to calculate free and bound concentrations of NADH from FLIM measurements is expected to hold in vivo. We next used this method to study NADH in mitochondria in oocytes. We applied Equations 2a,b to our FLIM data from oocytes and determined how the concentrations of free NADH, [ NADH f ] , and enzyme-bound NADH, [ NADH b ] , depended on oxygen level (Figure 1c, lower right). Interestingly, [ NADH f ] increased as oxygen fell below ~10 μM, while [ NADH b ] did not vary with oxygen level. We next explored the impact of metabolic inhibitors on mitochondrial NADH. We first inhibited LDH by adding 9 mM of oxamate to the AKSOM media. This led to a decrease of NADH intensity in the mitochondria (Figure 2a, upper) and significant changes in all FLIM parameters (Figure 2a, lower, p<0.001). We next inhibited complex I of the ETC using 5 μM of rotenone (in the presence of 9 mM of oxamate, to reduce cytoplasmic NADH signal for better mitochondrial segmentation). This resulted in a dramatic increase of NADH intensity in the mitochondria (Figure 2b, upper) and significant changes in NADH bound ratio and long lifetime (Figure 2b, lower, p<0.001). Then we inhibited ATP synthase with 5 μM of oligomycin (in the presence of 9 mM of oxamate), which, similar to rotenone, resulted in an increase of mitochondrial NADH intensity (Figure 2c, upper) and significant changes in all FLIM parameters (Figure 2c, lower, p<0.001). Finally, we subjected the oocytes to 5 μM of FCCP (in the presence of 9 mM of oxamate), which uncouples proton translocation from ATP synthesis, and observed a decrease of mitochondrial NADH intensity (Figure 2d, upper) and significant changes in FLIM parameters (Figure 2d, lower, p<0.001). Interestingly, the direction of change of FLIM parameters under FCCP is opposite to those under rotenone and oligomycin. For each of these conditions, we used Equations 2a,b to calculate the concentrations of free NADH, , (Figure 2f) from the measured intensity and FLIM parameters. While rotenone and oligomycin significantly increased . It remains unclear how to relate these changes of the free and bound concentrations of NADH to the activities of mitochondrial respiration.
Developing a coarse-grained NADH redox model to relate FLIM measurements of NADH to activities of mitochondrial metabolic pathways We next developed a mathematical model of NADH redox reactions to relate these quantitative FLIM measurements to activities of mitochondrial metabolic pathways. NADH is a central coenzyme that binds to enzymes and facilitates redox reactions by acting as an electron carrier. There are two categories of enzymes associated with NADH redox reactions, which together form a redox cycle: oxidases that oxidize NADH to NAD + and reductases that reduce NAD + to NADH. The major NADH oxidase in mitochondria is complex I of ETC for mammalian cells. There are many NADH reductases in mitochondria because NADH can be reduced through different pathways depending on the energy substrate. These pathways include the TCA cycle, fatty acid catabolism via beta oxidation, amino acid catabolism such as glutaminolysis and the malate-aspartate shuttle (Salway, 2017). A comprehensive NADH redox model will include all the oxidases and reductases. For generality, we consider N oxidases and M reductases.
For convenience, we introduced a reduced notation to describe models of the enzyme kinetics of these oxidases and reductases. We began by illustrating our reduced notation using reversible Michaelis-Menten kinetics as an example (Keleti, 1986;Miller and Alberty, 2002;Smith, 1992). The conventional, full notation for these kinetics (Figure 3a, left) explicitly displays all chemical species that are modeled in this reaction scheme -free NADH, NADH f , free enzyme, Ox i , free NAD + , NAD + f , and NADH bound to the enzyme -as well as the forward and reverse reaction ratesk −1 , k 1 , k −2 , and k 2 Our reduced notation for reversible Michaelis-Menten kinetics (Figure 3a, right) is an alternative way of representing the same mathematical model. In this reduced notation, only free NADH, free NAD + , and NADH bound to the enzyme are explicitly shown, while the free enzyme concentration is only represented as entering through the effective binding rates k b ox i and k ′ b ox i . The conventional, full notation, and the reduced notation are alternative ways of representing the same mathematical model, but the reduced notation is convenient to use in the derivation that follows (see Appendix 2).
We next introduced a generalized enzyme kinetics using our reduced notation (Figure 3b), which contains not only free NADH, free NAD + , and NADH bound to the enzyme, but also NAD + bound to the enzyme, and the reaction rates for oxidation and reduction of the bound coenzymes. In this reduced notation, all the binding and unbinding rates, and the reaction rates, can be functions of metabolite concentrations, protein concentrations, and other factors such as pH and membrane potential. As in the reversible Michaelis-Menten kinetics example, these rates can depend on the concentration of the free enzyme itself. This dependency on free enzyme concentration can be nonlinear, as could occur if the enzyme oligomerizes. Furthermore, the rates may depend on the concentration of free NADH, free NAD + , and the enzyme complexes. Thus, while the reduced notation for the generalized enzyme might appear to describe a first-order reaction, it can actually be used to represent reactions of any order, with arbitrary, nonlinear dependencies on the concentration of the enzyme itself, as well as arbitrary, nonlinear dependencies on other factors. In order to model the dynamics of enzymes described by such generalized kinetics, it is necessary to specify the functional form of all the rates, as well as specify mathematical models for all the variables that enter these rates (i.e., free enzyme concentration, membrane potential, pH, etc.) (Appendix 2). However, in what follows, we will derive results that hold true, irrespective of the functional form of the rates or the presence of additional, implicit variables. Thus, remarkably, these quantitative predictions are valid for enzyme kinetics of any order, with arbitrary nonlinearities in the rates.
To begin our derivation, we note that under this generalized enzyme kinetics (Figure 3b), the net flux through the ith oxidase at steady-state is: , and [ NADH f ] are the concentrations of the ith oxidase-bound NADH, NAD + , and free NADH, respectively. r + ox i and r − oxi are the forward and reverse oxidation rates. k b ox i and k u oxi are the binding and unbinding rates. The second equality in Equation 3 results from the steady-state condition, where the net binding and unbinding flux equals the net oxidation flux.
We next considered a redox cycle between NADH and NAD + with multiple oxidases and reductases. To account for all possible NADH redox pathways, we developed a detailed NADH redox model with N oxidases and M reductases described by the generalized enzyme kinetics (Figure 4a and Figure 4-figure supplement 1). In this model, NADH and NAD + can bind and unbind to each oxidase and reductase. Once bound, NADH can be reversibly oxidized to NAD + by the oxidases, and NAD + can be reversibly reduced to NADH by the reductases, forming a redox cycle. The functional dependencies of the binding and unbinding rates, and the reaction rates, can be different for each oxidase and reductase, and each of these rates can be nonlinear functions of free enzyme concentrations, NADH concentration, and other factors such as pH and membrane potential. Modeling the dynamics of this redox cycle requires specifying the precise number of oxidases and reductases, the functional forms of the rates, and mathematical models for all the variables the rates implicitly depend on. However, we will show that quantitative predictions regarding the interpretation of FLIM measurements can be made that generally hold, independent of these modeling choices.
FLIM cannot resolve the association of NADH with individual enzymes in cells, but rather, provides quantitative information on the global states of bound and free NADH. Thus, to facilitate comparison to FLIM experiments, we coarse-grained the detailed redox model by mapping all N oxidases into a single effective oxidase and all M reductases into a single effective reductase ( Figure 4b and Appendix 3). This coarse-graining is mathematically exact and involves no approximations or assumptions.
In the coarse-grained redox model, NADH can be bound to the effective oxidase, NADH · Ox , bound to the effective reductase, NADH · Re , or can be free, NADH f . Hence, the concentration of NADH bound to all enzymes is, ] , and the total concentration of NADH is, The kinetics of the effective oxidase and reductase are represented by the coarse-grained forward, r + ox , and reverse, r − ox , oxidation rates, and the forward, r + re , and reverse, r − re , reduction rates. The global flux through all the oxidases in the detailed redox model equals the global flux through the coarse-grained oxidase, which at steady-state is: where k b ox is the rate that free NADH binds the effective oxidase, k u ox is the rate that NADH unbinds the effective oxidase, and the last equality results because the coarse-grained redox loop is a linear pathway so the global oxidative flux must equal the global binding and unbinding flux at steady-state. The conservation of global flux explicitly relates the effective binding and unbinding rates and the reaction rates of the coarse-grained model to those of the detailed model (Appendix 3, Figure 4figure supplement 1). The binding and unbinding kinetics of NADH and NAD + to the effective oxidase and reductase are described by eight coarse-grained binding and unbinding rates (Figure 4b). The coarse-grained reaction rates and binding and unbinding rates can be arbitrary functions of metabolite concentrations, enzyme concentrations, and other factors (i.e., pH, membrane potential, etc.). These effective rates can even be functions of , and the concentration of other variables, and thus can include reactions of arbitrary order. Hence, this coarse-grained model is a generic model of NADH redox reactions. Fully specifying this model would require explicitly choosing the functional form of all the rates and incorporating additional equations to describe the dynamics of all the implicit variables that the rates depend on (Appendix 2). We next demonstrate that quantitative predictions regarding the interpretation of FLIM measurements of NADH can be made that are valid irrespective of the form of the rates or the presence of implicit variables.
Accurately predicting ETC flux from FLIM of NADH using the NADH redox model At steady-state, the model can be further coarse-grained, without approximation, to consider only free NADH, with a turnover rate of r ox , and free NAD + , with a turnover rate of r re (Figure 4c). Our key prediction is that the steady-state global oxidative flux of NADH is (Appendix 4): where and This prediction results from the steady-state assumption where the net binding and unbinding flux of NADH from the oxidase balances the net oxidative flux through the oxidase (Equation 4 and Appendix 4). The turnover rate of free NADH, r ox , is proportional to the difference between the NADH bound ratio β , that is, the ratio between bound and free NADH concentrations, and the equilibrium NADH bound ratio, βeq (i.e., what the bound ratio would be if the global oxidative flux is zero). βeq and the prefactor α are independent of the reaction rates of the oxidase and reductase and can be explicitly related to the binding and unbinding rates of the coarse-grained model (Appendix 4, Equation S43 and S45).
In mitochondria, the major NADH oxidation pathway is the ETC. Thus, Equations 5a-c predict that there is a direct connection between quantities that can be measured by FLIM of NADH in mitochondria (i.e., β and [ NADH f ] ) and the flux through the ETC (i.e., Jox ). Equations 5a-c suggest a procedure for using FLIM to infer flux through the ETC: if a condition can be found under which there is no net Free NADH, NADH f , and free NAD + , NAD + f , can bind and unbind with each oxidase and reductase. Once bound, NADH can be oxidized reversibly to NAD + by the oxidases, and NAD + can be reduced reversibly to NADH by the reductases, forming a redox cycle. Gray arrows represent the total fluxes through all oxidases and reductases of the redox cycle. (b) Coarse-grained NADH redox model. All oxidases and reductases are coarse-grained into a single effective oxidase and reductase, respectively. r + ox and r − ox are the coarsegrained forward and reverse oxidation rates of the oxidase; r + re and r − re are the coarse-grained forward and reverse reduction rates of the reductase.
re are the coarse-grained binding and unbinding rates of NADH and NAD + , respectively, to the oxidase and reductase. (c) At steady-state, all the kinetics of the model can be further coarse-grained into the turnover rate of free NADH, r ox , and the turnover rate of free NAD + , r re , characterizing the two branches of the cycle.
The online version of this article includes the following figure supplement(s) for figure 4: flux through the ETC, then βeq can be measured with FLIM. Once βeq is known, then subsequent FLIM measurements of β allows r ox , and hence Jox , to be inferred (up to a constant of proportionality α ) (Appendix 5).
Equations 5a-c are valid irrespective of the functional forms of the rate laws, which may have nonlinear dependencies on metabolite concentrations, enzyme concentrations, and other factors. While Equation 5a seems to imply first-order kinetics in [ NADH f ] , the rates can also be arbitrary functions of [ NADH f ] , so Equations 5a-c hold for kinetics of any order. Equations 5a-c are also applicable if the rates depend on additional variables that have their own dynamical equations (as long as the system is at steady-state): as an example, Appendix 7 shows that Equations 5a-c result when the N oxidases and M reductases are each described by reversible Michaelis-Menten kinetics, a model in which the rates depend on the concentration of free enzymes (which is a dynamical variable). More generally, if detailed biophysical models of the NADH oxidases are available, then parameters of these models can be explicitly mapped to the coarse-grained parameters of the NADH redox model. Appendix 7 and Appendix 7-table 1 contain mappings between the coarse-grained model and a number of previously proposed detailed biophysical models of NADH oxidation in the ETC (Beard, 2005;Korzeniewski and Zoladz, 2001;Hill, 1977;Jin and Bethke, 2002;Chang et al., 2011). However, since Equations 5a-c are valid for a broad set of models, they can be used for flux inference without the need to specify the functional form of the rates or the variables they depend on. This is because the rates are coarse-grained into two effective parameters α and βeq , which can be experimentally determined. This generality results from the steady-state assumption and the topology of the reactions resulting in the net binding and unbinding flux of NADH from the oxidase balancing the net oxidative flux. Thus, Equations 5a-c provide a general procedure to infer the ETC flux from FLIM measurements of NADH in mitochondria.
We applied this procedure to analyze our oxygen drop experiments ( Figure 1) by assuming that there was no net flux through the ETC at the lowest oxygen level achieved for each oocyte (implying that the measured value of β at that oxygen concentration corresponds to βeq for that oocyte). We also assumed that α and βeq do not change with oxygen levels, which is reasonable since, as noted above, they are independent of the reaction rates of the oxidase and reductase. The measured value of βeq allowed us to obtain a prediction for Jox as a function of oxygen concentration for the oocytes (Figure 5a). To test these predictions, we directly determined Jox as a function of oxygen concentration by measuring the OCR of the oocytes using a nanorespirometer (Lopes et al., 2005) (Materials and methods). The direct measurements of Jox from OCR quantitatively agree with the predictions of Jox from FLIM for all oxygen concentrations (Figure 5a), strongly arguing for the validity of the model and the inference procedure. This agreement supports the assumption that α and βeq are independent of oxygen levels.
So far, we have inferred the ETC flux up to a constant of proportionality α , allowing the relative changes of ETC flux to be inferred from FLIM of NADH. α cannot be determined by FLIM alone. If an absolute measurement of the ETC flux can be obtained at one condition, then α can be calibrated to predict absolute ETC fluxes for other conditions. OCR measurement provides a means to calibrate α (Appendix 5, Equation S48). We used oocytes cultured in AKSOM media at 50±2 μM oxygen as a reference state, which, from our OCR measurements yielded Jox = 56.6 ± 2 µM/s (SEM) and hence a constant of proportionality of α = 5.4 ± 0.2 s −1 . Using this value of α , we can predict absolute values of Jox under various perturbations assuming α remains a constant. We note that Jox is a flux density with units of concentration per second, an intensive quantity that does not depend on the mitochondrial volume. Multiplying Jox by the volume of mitochondria in an oocyte gives the total ETC flux, proportional to OCR, in that oocyte. In all subsequent discussions, ETC flux refers to flux density unless otherwise noted.
We next applied the inference procedure and a constant of α = 5.4 ± 0.2 s −1 to analyze the experiments of oxamate, FCCP, rotenone and oligomycin perturbations ( Figure 2). We dropped oxygen levels to determine βeq in the presence of oxamate ( Figure 5-figure supplement 1h) and applied Equations 5a-c to infer the impact of oxamate on Jox at 50 μM oxygen (i.e., control levels of oxygen). Surprisingly, while the addition of oxamate greatly impacts FLIM parameters, including a 29±2% (SEM) decrease in intensity and a 10±3% increase in bound ratio (Figure 2a), this procedure revealed that the predicted ETC flux with oxamate ( Jox = 55.2 ± 3.2 µM/s ) is the same as that without oxamate ( Jox = 55.4 ± 1.9 µM/s ) (Figure 5b; p=0.95), which was confirmed by direct measurements of oocytes' OCR that yielded Jox = 55.4 ± 1.5 µM/s and Jox = 54.9 ± 0.7 µM/s before and after the addition of oxamate, respectively (Figure 5b; p=0.85). We next analyzed the FCCP experiment. We obtained βeq by dropping oxygen in the presence of FCCP ( Figure 5-figure supplement 1h) and applied Equations 5a-c to infer the impact of FCCP on Jox at 50 μM oxygen. We predicted that FCCP increased the flux to Jox = 67.7 ± 1.5 µM/s , which was confirmed by the directly measured Jox = 74.0 ± 3.7 µM/s from OCR (Figure 5b; p=0.30). Following the same FLIM based inference procedures, we predicted that the addition of rotenone and oligomycin reduced the fluxes to Jox = 16.7 ± 1.4 µM/s and Jox = 15.7 ± 1.3 µM/s , respectively, which was again confirmed by corresponding direct measurements of OCR that yielded Jox = 11.1 ± 0.4 µM/s and Jox = 22.3 ± 0.6 µM/s (Figure 5b; p=0.31 and p=0.17). The quantitative agreement between predicted fluxes from FLIM and directly measured fluxes from OCR under a variety of conditions (i.e., varying oxygen tension, sodium oxamate, FCCP, rotenone, and oligomycin) demonstrates that Equations 5a-c can be successfully used to infer flux through the ETC in mouse oocytes. This agreement also supports the assumption that α is a constant across these different perturbations.
The work described above used the relation to predict the flux through the ETC from FLIM measurements. We next show that the model also predicts a relationship between r ox and the fluorescence lifetime of enzyme-bound NADH, τ l , in mitochondria. This provides a second means to use the model to infer r ox , and hence Jox , from FLIM of NADH. Specifically, we assumed that NADH bound to the oxidases have a different average lifetime, τox , than NADH bound to the reductases, τre , which is reasonable because NADH bound to different enzymes do exhibit different fluorescence lifetimes (Sharick et al., 2018). This assumption implies that the experimentally measured long lifetime of NADH in mitochondria, τ l , is a weighted sum of these two lifetimes, Using the coarse-grained NADH redox model at steady-state, Equation 6 leads to a non-trivial prediction that τ l is linearly related to 1/β (Appendix 5): where the slope A and offset B can be explicitly related to τox , τre , and the coarse-grained binding and unbinding rates. Such a linear relationship is indeed observed in individual oocytes subject to oxygen drops (Figure 6a and Figure 6-figure supplement 1), supporting the assumptions of the model. Combining Equations 7 and 5b leads to a predicted relationship between r ox and NADH long fluorescence lifetime (Appendix 5): where τeq is the equilibrium NADH long fluorescence lifetime, that is, the value of the long lifetime when the global oxidative flux is zero. This provides a second means to infer r ox from FLIM measurements: dropping oxygen and plotting the relationship between τ l and 1/β provides a means to measure A and B from Equation 7, while τeq can be obtained from the NADH long fluorescence lifetime obtained at the lowest oxygen level. Once A, B, and τeq are known, r ox can be inferred solely from NADH long fluorescence lifetime τ l , using Equation 8.
We next used the lifetime method (Equation 8) and the bound ratio method (Equation 5b) to separately infer r ox in oocytes subject to a wide variety of conditions (varying oxygen levels, with oxamate, FCCP, rotenone, and oligomycin). We obtained A, B, βeq , and τeq for these different conditions ( Figure 5-figure supplement 1 and Figure 6-figure supplement 1), and used the two different methods to provide two independent measures of r ox (assuming α is constant across all conditions). The predictions of r ox from these two methods quantitatively agree under all conditions (Figure 6b, p=0.73), which is a strong self-consistency check that further supports the use of the model to infer ETC flux from FLIM measurements of NADH. The NADH redox model enables accurate prediction of ETC flux in human tissue culture cells After thoroughly testing the NADH redox model and the inference procedure in mouse oocytes, we next investigated if it can be used in other cell types. We chose human tissue culture cells for this purpose, since they are widely used as model systems to study metabolic dysfunctions in human diseases including cancer (Vander Heiden et al., 2009) and neuropathology (Lin and Beal, 2006).
While mouse oocytes have a negligible level of NADPH compared to NADH (Bustamante et al., 2017), the concentrations of NADH and NADPH are similar in tissue culture cells (10-100 μM averaged over the whole cell) (Lu et al., 2018;Park et al., 2016;Blacker et al., 2014). Since NADPH and NADH have overlapping fluorescent spectra (Patterson et al., 2000), the presence of NADPH may complicate the interpretation of FLIM experiments. Thus, we investigated the impact of background fluorescence, such as from NADPH, on the flux inference procedure. If the background fluorescence does not change with the perturbations under study, then it can be treated as an additive offset that systematically makes the measured concentrations of free and bound NADH different from their actual values. In this case, a derivation in Appendix 5 demonstrates that the background fluorescence can be incorporated into the equilibrium bound ratio βeq and does not impact the flux inference procedures. In other words, if the modified βeq can be reliably determined, then the measured concentrations of free and bound fluorescent species can be used in place of the true values of NADH in Equations 5a-c to infer the ETC flux. An alternative possibility is that the background fluorescence does change with the perturbations under study, but in a manner that is proportional to the change in NADH. In this case, the background fluorescence can be incorporated into the equilibrium bound ratio α and, once   more, does not impact the flux inference procedures (Appendix 5). If the background fluorescence changes in some more complicated manner, then the inference procedure may no longer be valid. Thus, depending on the behavior of NADPH, it either might or might not interfere with the inference procedure: no impact if NADPH is either constant or proportional to changes in NADH, a possible impact otherwise. Therefore, the validity of the inference procedure in the presence of significant NADPH fluorescence must be established empirically. We next tested the inference procedures experimentally in hTERT-RPE1 (hTERT-immortalized retinal pigment epithelial cell line) tissue culture cells. We started by exploring the impact of metabolic perturbations on mitochondrial NAD(P)H: the combined signal from NADH and NADPH (which are indistinguishable) from mitochondria. We first cultured the cells in DMEM with 10 mM galactose (Materials and methods). We then inhibited complex I of the ETC by adding 8 μM of rotenone to the media. This resulted in a significant increase of mitochondrial NAD(P)H intensity (Figure 7a, upper). We segmented mitochondria using a machine learning-based algorithm from the intensity images of NAD(P)H, and fitted the fluorescence decay curves of mitochondrial NAD(P)H to obtain changes in FLIM parameters (Materials and methods). All FLIM parameters displayed significant changes (Figure 7a, lower, p<0.001, and Figure 7-figure supplement 1). We next uncoupled proton translocation from ATP synthesis by adding 3.5 μM CCCP to the media. This led to a decrease of NAD(P) H intensity in the mitochondria (Figure 7b, upper) and significant changes in NAD(P)H bound ratio and short lifetime, but in opposite directions as compared to rotenone perturbation (Figure 7b, lower, p<0.01, and Figure 7-figure supplement 1). Finally, we perturbed the nutrient conditions by culturing the cells in DMEM with 10 mM glucose. FLIM imaging revealed an increase of mitochondrial NAD(P)H intensity (Figure 7c, upper) and significant changes in all FLIM parameters as compared to the galactose condition (Figure 7c, lower, p<0.001, and Figure 7-figure supplement 1 The online version of this article includes the following source data and figure supplement(s) for figure 6: Source data 1. Excel spreadsheet of single-oocyte FLIM data used for Figure 6b.  Inference of the ETC flux from FLIM measurements requires a measurement of βeq . Since rotenone is known to drastically decrease the OCR of hTERT-RPE1 cells to near zero (MacVicar and Lane, 2014), we used the NAD(P)H bound ratio measured in the presence of rotenone as βeq . Different values of βeq were obtained for glucose and galactose conditions by adding 8 μM of rotenone to each condition (Figure 7-figure supplement 1). We next calculated the concentrations of free NAD(P)H, [ NAD(P)H f ] , from the FLIM parameters using Equation 2a.
[ NAD(P)H f ] displayed significant changes for all perturbations (Figure 7d). Using Equation 5b and assuming α is a constant, we calculated the NAD(P)H turnover rate, r ox , from the FLIM measurements and βeq . r ox changed significantly for all perturbations (Figure 7e). Multiplying r ox and [ NAD(P)H f ] , we obtained the predicted ETC flux, Jox , which increased under FCCP, decreased under glucose and reduced to zero under rotenone (Figure 7f).
To test the model predictions, we compared the predicted ETC flux with previous direct OCR measurements of the same cell type that we used, under the same conditions (MacVicar and Lane, 2014). Remarkably, the predicted changes in ETC fluxes are in quantitative agreement with the directly measured OCR across all conditions as estimated from Figure 1A of MacVicar and Lane, 2014: CCCP is predicted to increase the ETC flux by 14±3% (SEM), in agreement with the 18±21% increase from OCR measurement (p=0.80); Glucose is predicted to decrease ETC flux by 33±3%, in agreement with the 46±9% decrease from OCR measurements (p=0.30), shifting metabolism from oxidative phosphorylation to aerobic glycolysis. Since we used β from rotenone treatment as βeq , the predicted decrease in ETC flux after the addition of rotenone is 101±2%, which is in agreement with the 82±2% decrease from OCR measurement (p=0.28). This quantitative agreement between predicted ETC fluxes and measured OCR across all perturbations demonstrated the applicability of the NADH redox model and the flux inference procedures to tissue culture cells, even though they contain substantial levels of NADPH.
Homeostasis of ETC flux in mouse oocytes: perturbations of nutrient supply and energy demand impact NADH metabolic state but do not impact ETC flux Having established the validity of the NADH redox model and the associated flux inference procedures, we next applied it to study energy metabolism in mouse oocytes. We began by investigating the processes that determine the ETC flux in MII mouse oocytes. Mitochondrial-based energy metabolism can be viewed as primarily consisting of three coupled cycles: the NADH/NAD + redox cycle (which our NADH redox model describes), the proton pumping/dissipation cycle, and the ATP/ADP production/consumption cycle (Figure 8a). At the most upstream portion of this pathway, the reduction of NAD + to NADH is powered by a supply of nutrients, while at the most downstream portion, energy-demanding cellular processes hydrolyze ATP to ADP. To test whether nutrient supply and energy demand set ETC flux, we investigated the effect of perturbing these processes. To perturb supply, we first varied the concentration of pyruvate in the media from 181 μM (which is standard for AKSOM) to either 18.1 μM or 1.81 mM, and observed significant changes in NADH intensity and FLIM parameters (Figure 8b, left), demonstrating that the NADH metabolic state is altered. To perturb demand, we began by adding 10 μM nocodazole to the media, which disassembled the meiotic spindle, an energy user, and resulted in significant changes in NADH FLIM parameters (Figure 8b,  center). Similarly, the addition of 10 μM latrunculin A disassembled the actin cortex and also produced significant changes in NADH FLIM parameters (Figure 8b, right).
We next performed additional perturbations of nutrient supply, inhibiting the conversion of lactate to pyruvate by LDH (with 9 mM oxamate) and inhibiting the malate-aspartate shuttle (with 11 mM AOA). We performed additional perturbations of energy demand by inhibiting protein synthesis (with 1 mM cycloheximide) and ion homeostasis, by varying extracellular potassium concentrations from 0 mM to 15 mM, inhibiting the Na + /K + pump (with 2 mM ouabain), and adding an ionophore (10 μM gramicidin). All perturbations resulted in significant changes in NADH FLIM parameters (Figure 8figure supplement 1), showing that NADH metabolic state is generally impacted by varying nutrient supply and cellular energy demand. We next used the NADH redox model and the measured FLIM parameters to infer the concentration and effective turnover rate of free NADH for these perturbations. The free NADH concentrations, [ NADH f ] , and turnover rates, r ox , displayed large variations across the perturbations, ranging from 33.5±1 μM (SEM) to 56.0±2.8 μM and from 1.0±0.05 s -1 to 1.65±0.09 s -1 ,   respectively ( Figure 8c). Surprisingly, the changes in [ NADH f ] and r ox were highly anti-correlated such that the data points primarily fell within a region where the inferred ETC flux, Jox =rox [ NADH f ] , is a constant 55.5 μM·s -1 (Figure 8c, solid line, shaded region indicates 5% error). Indeed, ANOVA tests confirmed that perturbing nutrient supplies and cellular energy demand lead to no significant change in either the inferred ETC flux (Figure 8d, p=0.20) or directly measured OCR (Figure 8e, p=0.07). Thus, while nutrient supply and cellular energy demand strongly affect mitochondrial NADH redox metabolism, they do not impact ETC flux. In contrast, ETC flux is impacted by perturbing proton leak and ATP synthesis ( Figure 5). Taken together, this suggests that the ETC flux in mouse oocytes is set by the intrinsic properties of their mitochondria, which can adjust their NADH redox metabolism to maintain a constant flux when nutrient supplies and cellular energy demand are varied. The mechanistic basis of this homeostasis of ETC flux is unclear and will be an exciting topic for future research.
Subcellular spatial gradient of ETC flux in mouse oocytes: spatially inhomogeneous mitochondrial proton leak leads to a higher ETC flux in mitochondria closer to cell periphery Our results presented so far were performed by averaging together FLIM measurements from all mitochondria within an oocyte. However, FLIM data is acquired with optical resolution, enabling detailed subcellular measurements. To see if there are spatial variations in FLIM measurements within individual oocytes, we computed the mean NADH fluorescence decay time for each mitochondrial pixel. The mean NADH fluorescence decay time displays a clear spatial gradient, with higher values closer to the oocyte center (Figure 9a).
To quantify this gradient in more detail, we partitioned mouse oocytes into equally spaced concentric regions ( Figure 9b) and fitted the fluorescence decay curves from mitochondrial pixels within each region to obtain FLIM parameters as a function of distance from the oocyte center. NADH intensity, bound ratio, and long lifetime in mitochondria all display significant spatial gradient within oocytes ( Figure 9c). Next, using Equations 5a-c and βeq obtained at the lowest oxygen level, and confirming that βeq is uniform within the oocyte with complete inhibition of ETC using high concentration of rotenone (Figure 9-figure supplement 1), we predicted the ETC flux, Jox , as a function of distance from the oocyte's center. The ETC flux displayed a strong spatial gradient within oocytes, with a higher flux closer to the cell periphery (Figure 9d). Note that, as described above, Jox is actually a flux density with units of concentration per second. Thus, the measured flux gradient is not merely a reflection of variations in mitochondrial density, but instead indicates the existence of subcellular spatial heterogeneities in mitochondrial activities.
To investigate the origin of this flux gradient, we inhibited ATP synthase using 5 μM of oligomycin and repeated measurements of subcellular spatial variations in inferred fluxes. After inhibition, Jox decreased at all locations throughout the oocytes and displayed an even more dramatic flux gradient ( Figure 9e). If oligomycin completely blocks ATP synthase, then the remaining flux must be the result of proton leak. If it is further assumed that proton leak remains the same with and without oligomycin, then the flux due to ATP synthase in control oocytes can be determined by subtracting the flux after oligomycin inhibition (i.e., the proton leak) from the flux before inhibition. Performing this procedure throughout oocytes indicates that proton leak greatly increases in mitochondria near the periphery of oocytes, where ATP production decreases (Figure 9f). This implies that the subcellular gradient in ETC flux is primarily caused by a gradient in proton leak and that mitochondria near the periphery of oocytes are less active in ATP production than those in the middle of the oocyte.
We hypothesized that a gradient in proton leak would result in a gradient of mitochondrial membrane potential, with lower membrane potential closer to the cell periphery where proton leak is the greatest. To test this, we measured mitochondrial membrane potential using the membrane potential-sensitive dye TMRM, which preferentially accumulates in mitochondria with higher membrane potential (Al-Zubaidi et al., 2019). We observed a strong spatial gradient of the intensity of TMRM in mitochondria within oocytes, with dimmer mitochondria near the cell periphery (Figure 9g and h), indicating that mitochondria near the periphery of the oocyte have a lower membrane potential. This result is robust to locally normalizing TMRM intensity by mitochondrial mass using a membrane potential insensitive dye (Mitotracker Red FM), or using an alternative membrane potential-sensitive dye, JC-1 (Figure 9-figure supplement 2). The predicted flux of proton leak and mitochondrial TMRM intensity shows a strong negative correlation (Figure 9i), confirming our hypothesis. Figure 8. Homeostasis of ETC flux in mouse oocytes: perturbations of nutrient supply and energy demand impact NADH metabolic state but do not impact ETC flux. (a) The three coupled cycles of mitochondrial-based energy metabolism: the NADH/NAD + redox cycle, the proton pumping/ dissipation cycle, and the ATP/ADP production/consumption cycle. Nutrients supplied from the cytoplasm (blue) power the reduction of NAD + to NADH. Energy-demanding cellular processes in the cytoplasm (red) hydrolyze ATP to ADP. (b) Oocyte images (top) and change in NADH FLIM

Figure 8 continued on next page
Taken together, these results show that MII mouse oocytes contain subcellular spatial heterogeneities of mitochondrial metabolic activities. The observation that proton leak is responsible for the gradient of ETC flux suggests that the flux heterogeneity is a result of intrinsic mitochondrial heterogeneity. This is consistent with our conclusion from the homeostasis of ETC flux (Figure 8) that it is the intrinsic rates of mitochondrial respiration, not energy demand or supply, that controls the ETC flux. The causes and consequences of the subcellular spatial variation in mitochondrial activity remain unclear and are an exciting topic for future research.

Discussion The NADH redox model is a general model to relate FLIM measurements of NADH to ETC fluxes
Despite extensive studies and applications of FLIM in metabolic research (Bird et al., 2005;Skala et al., 2007;Heikal, 2010;Sharick et al., 2018;Sanchez et al., 2018;Liu et al., 2018;Sanchez et al., 2019;Ma et al., 2019), it remains a challenge to relate FLIM measurements to the activities of the underlying metabolic pathways in cells. We overcame this challenge by developing a coarsegrained NADH redox model that leads to quantitative predictions for the relationship between FLIM measurements and the flux through the ETC. The model was constructed by explicitly coarse-graining a detailed NADH redox model with an arbitrary number of oxidases and reductases that represent all the possible enzymes involved in NADH redox reactions. The reactions in the detailed NADH redox model can be of arbitrary order and depend on implicit variables (i.e., free enzyme concentration, membrane potential, pH, etc.), which obey their own dynamical equations. The dynamics of the redox model will, of course, depend on the precise number of oxidases and reductases, the functional forms of the rates, and specific mathematical models for all the variables the rates implicitly depend on. However, the quantitative predictions relating FLIM measurements and ETC flux are independent of these modeling choices. Coarse-graining the detailed NADH redox model reduces all oxidases to an effective oxidase and all reductases to an effective reductase. The kinetic rates of the coarsegrained model can be related to those of the detailed model by keeping the global fluxes through the oxidases and the reductases the same in both models. The coarse-grained model predicts that the flux through the ETC is a product of the turnover rate and the concentration of free NADH (Equation 5a). The turnover rate is proportional to the difference between the nonequilibrium and the equilibrium NADH bound ratio (Equation 5b), which are measurable by FLIM of NADH (Equation 5c). Thus, this model provides a generic framework to relate FLIM measurements of NADH to the flux through the ETC in mitochondria.
The central assumption required for the validity of Equations 5a-c is that the redox reactions, and binding and unbinding processes, can be approximated as being at steady-state (i.e., undergoing only quasistatic changes over perturbations or development). At steady-state, the net binding and unbinding flux balances the oxidative flux of NADH. Therefore, the measurement of binding and unbinding state of NADH from FLIM allows the inference of the ETC flux, irrespective of the detailed behaviors of the oxidative reactions. parameters relative to control (bottom) for changing pyruvate concentration (left), addition of 10 μM nocodazole (center) and addition of 10 μM latrunculin A (right). Student's t-test were performed for the change of FLIM parameters (*p<0.05, **p<0.01, ***p<0.001). The spindle disassembles after addition of 10 μM nocodazole (top, center) and the actin cortex disassembles after addition of 10 μM latrunculin A (top, right). (c) NADH turnover rate ( r ox ) and NADH free concentrations ( ) inferred from FLIM measurements under a variety of perturbations of nutrient supply and energy demand. Error bars are standard error of the mean (s.e.m) across oocytes. The black line corresponds to r ox and values with an inferred flux of Jox = 55.5 µM · s −1 , and the gray shaded region corresponds to a variation of ±5% around that value. (d) The inferred ETC flux and (e) measured OCR show no change across different perturbations of nutrient supply and energy demand (ANOVA, p=0.20 and p=0.07, respectively). ETC, electron transport chain; FLIM, fluorescence lifetime imaging microscopy; OCR, oxygen consumption rate.
The online version of this article includes the following source data and figure supplement(s) for figure 8: Source data 1. Excel spreadsheet of single-oocyte FLIM data and batch OCR data used for    Remarkably, all the binding and unbinding rates of the NADH redox model are coarse-grained into two effective parameters: α and βeq , which can be experimentally measured. We determined the value of α from an OCR measurement (Appendix 5, Equation S48), and we determined the value of βeq from FLIM of NADH at low oxygen levels or from rotenone perturbation (Appendix 5, Figure 5-figure supplement  1h). In MII mouse oocytes, α do not significantly vary in response to oxygen or drug and nutrient perturbations. This is demonstrated by the agreement between the predicted ETC flux and the measured OCR with a constant α of 5.4±0.2 s -1 across a variety of conditions ( Figure 5). α is predicted to depend only on the coarse-grained unbinding rates of NADH from the enzymes (Equation S43), so the observed constancy of α implies that the perturbations in this study primarily impacted the reduction/oxidation reaction rates (and not the unbinding rates). In other scenarios, such as when the concentrations of enzymes change, the coarse-grained unbinding rates might change, so α might not be a constant. In contrast, βeq does vary with drug and nutrient perturbations, but not with oxygen level, allowing βeq to be obtained at the lowest oxygen level for different drug and nutrient conditions ( Figure 5-figure supplement 1h and Figure 8-figure  supplement 1). Using these two parameters, we inferred the effective turnover rate of free NADH, r ox , from FLIM measurements of NADH. By multiplying this turnover rate with the concentration of free NADH, [ NADH f ] (also obtained from FLIM measurements using Equation 2a), we inferred the ETC flux from Equation 5a. Thus, all the complex behaviors of the binding and unbinding and reaction rates are captured by the variations in FLIM parameters of NADH, and our coarse-grained model provides a generic way to interpret these variations.
While we found that βeq is smaller than β in mouse oocytes, this does not generically have to be true. Thus, if a perturbation is observed to decrease the NADH bound ratio β, it does not necessarily imply a decrease of the ETC flux. Similarly, a decrease of NADH long lifetime is not necessarily associated with an increase of the ETC flux. Therefore, measurements of α and βeq are required to use Equations 5a-c to infer ETC flux from FLIM measurements of NADH.

The underlying assumptions and limitations of the NADH redox model
In this section, we clarify the underlying assumptions and limitations of the model to facilitate the accurate interpretation of FLIM measurements of NADH in different biological contexts.
To use the coarse-grained NADH model, segmentation needs to be performed to separate the mitochondrial NADH signal from the cytoplasmic NADH signal, because they encode different metabolic fluxes. In mouse oocytes, the segmentation can be reliably performed based on NADH images due to the higher NADH intensity in mitochondria than cytoplasm. Mitochondrial movements are also slow in MII oocytes (Video 1); hence, long exposure times can be used to obtain high contrast NADH images. For cells where NADH contrast is low, such as in yeast cells (Papagiannakis et al., 2017;Shaw and Nunnari, 2002), MitoTracker dye (Appendix 1 Figure 1-figure supplement 1) or mitochondrial associated fluorescent proteins (Westermann and Neupert, 2000) will likely be needed for reliable segmentation of mitochondria.
One of the most important assumptions that enables the coarse-grained model to be used to predict fluxes is that the NADH redox cycle can be well approximated as being at steady-state, that is, the rate of change of NADH concentrations is much slower than the kinetic rates, including the binding/unbinding rates and the reaction rates. This is true for mouse oocytes, where the NADH intensity does not significantly change over the course of hours. This assumption also holds for slow processes such as the cell cycle (Papagiannakis et al., 2017), which occurs on the timescale of hours from the oocyte center (n=67). (d) Predicted ETC flux from FLIM of NADH as a function of distance from the oocyte center (n=67). (e) ETC flux gradient is enhanced by 5 µM oligomycin (n=37), suggesting the flux gradient is determined by proton leak. CT is AKSOM with oxamate (n=32). 9 mM oxamate is present in oligomycin condition to reduce cytoplasmic NADH signal for better mitochondrial segmentation. (f) Opposing flux gradients of proton leak and ATP production, where proton leak (ATP production) is maximal (minimal) at the cell periphery. (g) Heatmap of the TMRM intensity in mitochondria, which increases with mitochondrial membrane potential, exhibits a subcellular spatial gradient within oocytes. (h) Mitochondrial TMRM intensity as a function of distance from the oocyte center (n=16). (i) Predicted flux of proton leak correlates negatively with mitochondrial membrane potential as measured by mitochondrial TMRM intensity. Scale bar, 20 μm. Error bars represent standard error of the mean (s.e.m) across different oocytes. ETC, electron transport chain; FLIM, fluorescence lifetime imaging microscopy.
The online version of this article includes the following figure supplement(s) for figure 9:   compared to timescales of seconds for the kinetic rates. This claim is supported by the success of the model on human tissue culture cells. The steady-state approximation could fail for rapid dynamics of NADH, such as the transient overshoot of NADH in neurons induced by acute external stimulus (Díaz-García et al., 2021), but this needs to be tested experimentally.
While NADH and NADPH share the same fluorescence spectrum, NADH concentration is 40 times greater than the concentration of NADPH for the whole mouse oocytes and presumably even higher for mitochondria (Bustamante et al., 2017). NADPH concentration can be comparable to that of NADH for other cell types such as tissue culture cells (Park et al., 2016). However, we have shown that the presence of NADPH signal and other background fluorescence signals only affect the equilibrium bound ratio βeq or the prefactor α , and hence does not affect the flux inference procedure if βeq can be reliably determined and α remains a constant (Appendix 5). This was validated in tissue culture cells by comparing predicted ETC flux (Figure 7) with previous OCR measurements ( MacVicar and Lane, 2014).
Finally, when relating NADH FLIM measurements to the ETC flux, we did not explicitly consider the contribution to the flux through FADH 2 . This is a valid approximation when the FADH 2 oxidative flux is much smaller than the NADH oxidative flux, as is often the case since pyruvate dehydrogenase plus the TCA cycle yields four NADH molecules but only one FADH 2 molecule per cycle. Alternatively, if the FADH 2 flux is proportional to the NADH flux, then a rescaled value of α can be used in Equation 5b to effectively account for both fluxes. The proportionality of FADH 2 flux and NADH flux is expected when NADH and FADH 2 are produced from the same redox cycle with fixed stoichiometry, such as the pyruvate dehydrogenase and TCA cycle. This proportionality will break down if significant amounts of NADH and FADH 2 are produced in independent cycles where the stoichiometry varies, for example, when the glycerol phosphate shuttle acts as a reductase in mitochondria for FADH 2 but not for NADH.
Given these underlying assumptions, the model needs to be tested before being applied to other biological systems. The present study provides an example for such tests in mouse oocytes and human tissue culture cells by comparing the predicted ETC flux from FLIM with direct measurements of OCR across a wide range of perturbations.

Towards spatiotemporal regulations of metabolic fluxes in cells
Cells transduce energy from nutrients to power various cellular processes. The ETC flux represents the total rate of energy conversion by mitochondria. Despite the detailed knowledge of the biochemistry of mitochondrial metabolism, it is still unclear what cellular processes determine ETC flux or how cells partition energetic fluxes to different cellular processes, including biosynthesis, ion pumping, and cytoskeleton assemblies. Energetic costs of specific cellular processes have been estimated from theoretical calculations (Stouthamer, 1973) or through inhibition experiments (Mookerjee et al., 2017). The latter typically involves measurements of the change of metabolic fluxes, such as OCR, upon inhibition of specific cellular processes, and interpreting this change as the energetic cost of the inhibited process. This interpretation is valid if metabolic flux is determined by the energy demand of different cellular processes in an additive manner. This assumption has not been thoroughly tested. Using the NADH redox model, we discovered a homeostasis of ETC flux in mouse oocytes where perturbing energy demand and supply do not impact ETC flux despite significantly changing NADH metabolic state. On the other hand, perturbing ATP synthesis and proton leak greatly impacted the ETC flux. From these results, we concluded that it is the intrinsic rates of mitochondrial respiration, rather than energy supply or demand, that controls the ETC flux in mouse oocytes. While NADH metabolic state significantly changed in response to perturbing energy demand and supply, indicating cell metabolism was indeed impacted, it is unclear if these perturbations also influenced ATP, ADP, or AMP levels. Future work, including direct measurements of ATP, ADP, and AMP levels, will be required to uncover the mechanism of flux homeostasis. More broadly, our work demonstrates that it is a prerequisite to understand the regulation of ETC fluxes in order to correctly interpret the changes of ETC flux upon inhibiting subcellular processes.
The mechanism of the homeostasis of ETC flux is unclear. One possibility is the presence of flux buffering pathways, where the change of ATP fluxes induced by process inhibition is offset by the opposing change of fluxes through the buffering pathways. Enzymes such as adenylate kinase are known to buffer concentrations of adenine nucleotide (De la Fuente et al., 2014), but it is unclear if they also buffer fluxes. Another possibility is a global coupling of cellular processes, where the change of ATP consumption by one process is offset by the change of others. Changes in proton leak could also compensate for changes in ATP production. Additional work will be required to distinguish between these (and other) possibilities.
FLIM data is obtained with optical resolution, enabling subcellular measurements of NADH metabolic state. Interpreting these measurements using the NADH redox model enables inference of metabolic fluxes with subcellular resolution. Using this method, we discovered a subcellular spatial gradient of ETC flux in mouse oocytes, where the ETC flux is higher in mitochondria closer to the cell periphery. We found that this flux gradient is primarily a result of a spatially heterogeneous mitochondrial proton leak. It will be an exciting aim for future research to uncover the causes and consequences of the subcellular spatial variation in mitochondrial activity.

Culturing of mouse oocytes
Frozen MII mouse oocytes (Strain B6C3F1) were purchased from EmbryoTech. Oocytes were thawed and cultured in droplets of AKSOM media purchased from MilliporeSigma in plastic petri dish. Mineral oil from VitroLife was applied to cover the droplets to prevent evaporation of the media. Oocytes were then equilibrated in an incubator at 37°C, with 5% CO 2 and air saturated oxygen before imaging. For imaging, oocytes were transferred to a 2-μl media droplet in a 35-mm glass bottom FluoroDish from WPI covered with 400-500 μl of oil. The glass bottom dish was placed in an Ibidi chamber with temperature and gas control during imaging. Temperature was maintained at 37°C via heated chamber and objective heater. CO 2 was maintained at 5% using gas tanks from Airgas.

Cell lines
The hTERT-RPE1 cell line is an established wild-type cell line received from the Cheeseman lab that has been validated based on behavior and properties. The hTERT-RPE1 cell line was maintained and tested for mycoplasma contamination in the Needleman lab on a regular basis (Southern Biotech).

Culturing of hTERT-RPE1 cells
Cell lines were maintained at 37°C and 5% CO 2 . Cells were grown in Dulbecco's modified Eagle's medium (DMEM) (11966025, Gibco) supplemented with 10% fetal bovine serum (FBS), 0.5 mM sodium pyruvate, 5 mM HEPES, 1% penicillin and streptomycin, and either 10 mM glucose or 10 mM galactose. Cells were passaged in glucose or galactose at least three times before imaging. Cells were plated on 35 mm glass bottom FluoroDishes from WPI for imaging. Right before imaging, the media was replaced with 1 ml of phenol red-free DMEM (A1443001, Gibco) supplemented with 0.5 mM sodium pyruvate, 4 mM L-glutamine, 10 mM HEPES, and either 10 mM glucose or 10 mM galactose.

FLIM measurements
Our FLIM system consists of a two-photon confocal microscope with a 40× 1.25 NA water immersion Nikon objective, Becker and Hickle Time Correlated Single Photon Counting (TCSPC) acquisition system and a pulsed MaiTai DeepSee Ti:Sapphire laser from Spectra-Physics. NADH autofluorescence was obtained at 750 nm excitation wavelength with a 460/50 nm emission filter. Laser power at the objective was maintained at 3 mW. The scanning area was 512 by 512 pixels with a pixel size of 420 nm. Acquisition time was 30 s per frame. Oocytes were imaged with optical sectioning across their equators. A histogram of NADH fluorescence decay times was obtained at each pixel of the image.

Oxygen measurements
Oxygen level was measured in the Ibidi chamber with an electrode-based oxygen sensor (GasLab).
Since the oil layer covering the media droplet was very thin, the oxygen level in the droplet was assumed to be in instant equilibration with the chamber.

Image and FLIM data analysis
To separate mitochondrial NADH signal from cytoplasmic signal, we performed machine learningbased segmentation algorithms on NADH intensity images. We used the freeware Ilastik (Berg et al., 2019), which implements a supervised learning algorithm for pixel classification. The classifiers were trained to separate mitochondrial pixels from cytoplasmic pixels with a greater than 80% accuracy, as tested by MitoTracker Red FM (Appendix 1, Figure 1-figure supplement 1). We grouped photons from all mitochondrial pixels to obtain a histogram of NADH decay times for each oocyte and for each image of tissue culture cells. To extract the FLIM parameters of NADH bound fraction f, long lifetime τ l and short lifetime τs , we fitted the histogram with G = IRF * ( , where * indicates a convolution, and IRF is the instrument response function of the FLIM system, measured using a urea crystal. is the two-exponential model for the NADH fluorescence decay. C 1 is the amplitude of the decay and C 2 is the background noise. The fitting was performed with a custom MATLAB code using a Levenberg-Marquardt algorithm (Yoo, 2018). To obtain the intensity, I , of mitochondrial NADH, we first measured the average number of photons per mitochondrial pixel, and divided it by the pixel area, 0.185 μm 2 , and pixel scanning time 4.09 μs. The flux of ETC is inferred using Equations 5a-c for each oocyte and for tissue culture cells in a single image. Heatmaps of mean NADH fluorescence decay times were obtained by computing NADH fluorescence decay time of each mitochondrial pixel and averaging over neighboring mitochondrial pixels weighted by a Gaussian kernel with a standard deviation of 20 pixels. All FLIM measurements were taken from distinct individual oocytes and distinct images of tissue culture cells. Error bars in all figures of FLIM represent standard error of the mean across different individual oocytes or across different images for tissue culture cells. Number of oocytes is reported with n. Number of images for tissue culture cells is reported with N.

Error analysis
FLIM curves were independently fit for each individual oocyte. The reported error bars in this manuscript are standard errors of the mean (SEMs) across these measurements, which depends on the level of variation (the standard deviation) between the oocytes. Two sources of variation in FLIM measurements across the oocytes are: (1) true biological variations between oocytes and (2) fitting errors in the FLIM analysis. To estimate the error of fitting, we performed bootstrapping with randomly drawn points with substitution from each fluorescence decay curve for 53 oocytes. There are ~66,000 photons per oocyte, from which we generated 10 bootstrapped decay curves per oocyte to estimate the fitting error. The fitting error is computed as the variance and covariance of the fitted parameters across bootstrapped decay curves and averaged over 53 oocytes. At high oxygen level in the AKSOM condition, the bootstrapping yields a variance of 2.2×10 -4 , 4.6×10 -3 ns 2 , and 6.0×10 -4 ns 2 for bound fraction, long lifetime, and short lifetime, respectively. The cell-to-cell variances obtained from a single fit per oocyte are 4.4×10 -4 , 9.5×10 -3 ns 2 , and 1.6×10 -3 ns 2 for bound fraction, long lifetime, and short lifetime, respectively. Hence the bootstrapping error accounts for 50%, 49%, and 40% of the cell-to-cell variance in bound fraction, long lifetime, and short lifetime, respectively. The bootstrapping yields a covariance of -1.0×10 -3 ns between bound fraction and long lifetime, which only accounts for ~20% of the covariance between these two variables during oxygen drop experiment. The inferred mean flux for oocytes at high oxygen levels in AKSOM is ⟨Jox⟩ = 56.6 µM · s −1 . Propagating the error of fitting in all parameters from the bootstrapping analysis to the inferred flux gives a standard error of the mean in Jox of 1.1 μM·s -1 . The standard error of the mean in Jox obtained from a single fit per oocyte was 2.0 μM·s -1 . Thus, fitting errors account for ∼50% of the standard error of the mean in Jox .

Metabolic and demand perturbations
Oxygen drop experiments for oocytes were performed by mixing nitrogen-balanced 5% O 2 gas with 0% O 2 gas at different ratios to create a continuous oxygen drop profile. CO 2 was maintained at 5%. Oocytes were imaged for 10 min at 5% O 2 , 30 min during the continuous drop from 5% O 2 to approximately 0% O 2 , and 20 min after quickly returning to 5% O 2 . Oxygen levels were simultaneously monitored with an electrode-based oxygen sensor in the Ibidi chamber. 5% O 2 corresponds to ~50 μM of oxygen concentration in the culturing media. All the drug perturbations for oocytes were performed by equilibrating oocytes in the AKSOM media containing the corresponding drug for 15-30 min before the oxygen drop experiments. Pyruvate and potassium perturbations were performed by making KSOM media following Cold Spring Harbor Laboratory protocols with varying concentrations of sodium pyruvate and potassium, respectively. For oligomycin, FCCP, rotenone and pyruvate perturbations, 9 mM of sodium oxamate was also added to the media to suppress cytoplasmic NADH signal for better mitochondrial segmentation. The addition of the oxamate does not change the ETC flux of the mitochondria (Figure 5b).
For hTERT-RPE1 cells, drug perturbations were performed by replacing the media with drug-containing media through pipetting. Cells were imaged for 20-30 min immediately after drug perturbations.
All drugs were purchased from Sigma-Aldrich. Temperature was maintained at 37°C. CO 2 was maintained at 5%.

Oxygen consumption rate measurement
The OCR of the oocytes was measured using the nanorespirometer from Unisense (Lopes et al., 2005). A batch of 10-15 oocytes was placed at the bottom of a glass capillary with a diameter of 0.68 mm and a height of 3 mm. The capillary well is filled with AKSOM media or drug-containing media for metabolic perturbations. After an equilibration time of ~2 hr, a steady-state linear oxygen gradient is established in the capillary well due to the balance of oocyte respiration and oxygen diffusion. A motor-controlled electrode-based oxygen sensor (Unisense) is used to measure the oxygen gradient. The OCR is calculated as the product of the oxygen gradient, diffusivity of oxygen in the media, taken to be 3.37×10 -5 cm 2 /s, and the cross-sectional area of the capillary well, which was 0.36 mm 2 . The entire system was enclosed in a custom-built chamber with temperature and gas control. Temperature was maintained at 37°C. Oxygen level was continuously varied during oxygen drop experiments by slowly mixing 20% O 2 with 0% O 2 from gas tanks, and maintained at the air saturation level for drug and pyruvate perturbations. OCR was measured on a group of 10-15 oocytes at a time. Single-oocyte OCR was obtained by dividing the measured OCR by the number of oocytes in the group. Error bars in all figures of OCR represent standard error of the mean across different groups of oocytes normalized by the number of oocytes in each group. Number of oocytes is reported with n. Number of groups is reported with N.

Statistical analysis
For the comparison between inferred ETC flux and measured ETC flux of the oocytes, two-sample t-test was performed on the vectors of inferred single-cell ETC flux (with n elements, where n is the number of oocytes) and the batch OCR measurements (with N elements, where N is the number of batch groups). For the comparison between inferred ETC flux and measured ETC flux of the tissue culture cells, twosample t-test was performed on the vectors of inferred relative change of ETC flux (with n elements, where n is the number of images) and the relative change of OCR estimated from Figure 1A of MacVicar and Lane, 2014 (with N elements, where N is the estimated number of OCR data points).

Mitochondrial membrane potential measurement
The spatial distribution of mitochondrial membrane potential within oocytes was measured with a potential-sensitive dye TMRM (Sigma-Aldrich). Oocytes were cultured in AKSOM with 100 nM TMRM for 30 min before imaging. TMRM signal was obtained at 830 nm excitation wavelength with 560/40 nm emission filter. Mitochondrial TMRM intensity in different regions of the oocyte was computed by dividing the total number of photons from that region by the number of pixels in the same region. Heatmaps of mitochondrial TMRM intensity were obtained by computing photon counts for each mitochondrial pixel and averaging over neighboring mitochondrial pixels weighted by a Gaussian kernel with a standard deviation of 20 pixels. To normalize TMRM intensity by mitochondrial mass, we cultured oocytes in AKSOM with 100 nM MitoTracker Red FM and 25 nM TMRM for 30 min before imaging. We also cultured oocytes in AKSOM with 1 μg/ml JC-1 dye for 3 hr before imaging.
Mitochondrial membrane potential of hTERT-RPE1 cells was measured with TMRM. The cells were cultured in DMEM with 100 nM TMRM for 15-30 min before imaging. To measure membrane potential under drug perturbations, the original media was pipetted out and replaced with media containing both 100 nM TMRM and the drug. The cells were imaged for 20-30 min immediately after drug perturbations. TMRM intensity ratio was obtained by normalizing the mitochondrial TMRM intensity by the cytoplasmic TMRM intensity.

Data availability
All data generated or analysed during this study are included in the manuscript and supporting file; Source Data files have been provided for Figures 2, Figure

Appendix 1 Segmentation of mitochondria and calculation of NADH concentrations Segmentation of mitochondria
We used Ilastik, a machine learning-based software for image analysis, to classify pixels in the NADH intensity images containing mitochondria (Berg et al., 2019). For each experiment, we generated a time-lapse movie of NADH (Video 1). We used a few images in the movie as the training data set to train the software to classify mitochondrial pixels by manually selecting clustered high brightness pixels. Other pixels are classified as either cytoplasm or background. We then applied the trained pixel classifier to generate a mitochondrial probability map for each image in the entire movie, with each pixel assigned a probability between 0 and 1 to be mitochondrial pixel. Pixels with a probability higher than 0.7 were considered to be mitochondrial pixels.
To test the accuracy of this segmentation algorithm, we immersed the oocytes in AKSOM media containing MitoTracker Red FM, a dye that specifically labels mitochondria. Pixels with intensity above 60 percentile in the MitoTracker image were considered to be mitochondrial pixels. We imaged NADH and MitoTracker for the same oocyte and compared the resulting distribution of mitochondria (Figure 1-figure supplement 1). We defined the accuracy of the NADH-based segmentation as the fraction of photons originating from true mitochondrial pixels. The accuracy of the segmentation is 80.6±1% (SEM) for the control condition as averaged over seven oocytes. We repeated the analysis for oxamate, oligomycin, FCCP and rotenone perturbations, and obtained an accuracy of segmentation of 78.6±1.4%, 84.1±1.6%, 83.7±0.5%, and 81.7±2%, respectively, similar to the control condition.

Converting NADH intensity to NADH concentrations
Since the molecular brightness of NADH depends on the fluorescence lifetime of NADH, which changes drastically upon binding enzymes, the NADH concentration is not linearly proportional to NADH intensity. FLIM provides an accurate way of measuring NADH concentrations by simultaneously measuring fluorescence intensity and lifetime. We now derive the NADH intensityconcentration relation from the FLIM measurements. Assuming molecular brightness is proportional to the fluorescence lifetime, and therefore that free and bound NADH have different contributions to the measured intensity, we have where I is the intensity of NADH and cs is a calibration factor that depends on the laser power. From Equation S1, we obtained the concentrations of free and bound NADH: where f is the fraction of bound NADH.
To get the calibration factor cs , we titrated NADH in AKSOM solutions and fitted the calibration curve using: where τ sol is the lifetime of NADH in solution. τ sol was directly measured by FLIM (Figure 1-figure  supplement 2), allowing us to obtain cs from the fit (Figure 1-figure supplement 2).

FLIM can be used to accurately measure concentrations of bound and free NADH in vitro
To test if absolute concentrations of free and bound NADH can be accurately measured from FLIM of NADH, we prepared solutions with known total concentration of NADH, and titrated the concentration of LDH, an enzyme to which NADH can bind. We prepared the solutions with 50 mM TRIS buffer, 150 mM NaCl at pH 7.6 and 37°C. We added a total concentration of 50 μM NADH to the solution and titrated LDH concentrations at 0, 1.4, and 3.5 μM. We first performed single exponential fitting of the NADH decay curve at 0 μM LDH, where all NADH are free, to obtain the NADH short lifetime (Figure 1-figure supplement 3d). From the NADH intensity (Figure 1-figure  supplement 3a), we obtained the calibration factor cs using Equation S4 with [ NADH sol ] = 50µM . We then fixed the short lifetime and performed two-exponential fitting of the NADH decay curve at LDH concentrations of 1.4 μM and 3.5 μM to obtain the bound ratio (Figure 1-figure supplement 3b) and long lifetime (Figure 1-figure supplement 3c). As expected, NADH bound ratio increases with LDH concentrations, as there is more enzyme for NADH to bind. Finally, we calculated free NADH concentration [ NADH f ] and bound NADH concentration [ NADH b ] using Equations S2 and S3 from the FLIM parameters. Remarkably, the free and bound concentrations of NADH both change with LDH concentrations but the total concentration remains at 50 μM (Figure 1-figure supplement  3e). This result shows that Equations S2-S3 can be used to accurately measure the concentrations of free and bound NADH from FLIM measurements of NADH.
effective binding and unbinding rates of NADH to the oxidase, k ′ b ox i and k ′ u ox i are the effective binding and unbinding rates of NAD + to the oxidase, and r + ox i and r − oxi are the forward and reverse oxidation rates. These rates can be arbitrary functions of implicit variables, such as the concentration of free oxidase, [ Ox i ] , the mitochondrial membrane potential, ∆G H , pH , and other factors: Thus, while Equations S13-S16 superficially appear to be linear and first order, they can actually refer to nonlinear reactions of any order because the rate can depend on [ NADH f ] and other variables (Equations S17a-f).
The implicit variables that these rates depend on can each be governed by their own dynamics that are arbitrary functions of other variables: Describing the dynamics of the enzyme requires specifying the implicit variables that the rates depend on, the functional form of these dependencies, and the additional equations for the dynamics of the implicit variables (Equations S18). However, we will show that the predicted relationship between FLIM measurements of NADH and fluxes does not depend on these modeling choices. Thus, the reduced notation is convenient for deriving these relations for a broad class of models. which leads to We require the global forward and reverse reaction flux through the effective oxidase and reductase to be equal to the sum of the reaction fluxes through all of the individual oxidases and reductases: which leads to By applying the same procedure to NAD + , we can obtain the effective reduction rates r + re , r − ox and the effective binding and unbinding rates of NAD + : k We omit the derivation here because these rates are not needed to infer ETC flux. We hence explicitly related the kinetic rates of the coarse-grained model (Figure 4b) to those of the detailed model (Figure 4-figure supplement  1). We note that under the generalized enzyme kinetics (Figure 3b; Equations S13-S16), all kinetic rates are considered to be general functions of enzyme concentrations, metabolite concentrations, and other factors, and thus, all of the rates in the coarse-grained model can also depend on all of those factors. These implicit variables can obey their own dynamical equations (Equation S18). The coarse-graining presented here is mathematically exact and independent of both the functional forms of these rates and the functional form of the dynamic equations of the implicit variables.
where we defined the NADH bound ratio and its equilibrium counterpart as: (S47) , provides a method to obtain r ox . We measure the NADH bound ratio, β , using β = f /(1 − f) , where f is the NADH bound fraction obtained by fitting the fluorescence decay curve of NADH (see Materials and methods). We obtain the equilibrium bound ratio, βeq , by dropping the oxygen level to the lowest achievable value with our setup: [ = 0.26 ± 0.04 µM , and assuming βeq does not change with oxygen levels. Note that βeq does change with drug perturbations, and therefore needs to be separately determined for each condition ( Figure 5-figure supplement  1h). We obtain α using direct measurement of Jox from OCR measurements: where Vm = 9.5 × 10 4 µm 3 is the average volume of mitochondria per oocyte approximated from the area fraction of mitochondria based on the segmentation, where the mitochondrial area fraction is estimated at 46% and oocyte volume at 2×10 5 µm 3 . Using OCR = 2.68 ± 0.06 fmol/s per oocyte in the control condition (AKSOM media at 50 μM oxygen level), we get α = 5.4 ± 0.2 s −1 . α is approximated as a constant that does not vary with perturbations, hence α calibrated at one condition can be used for all other conditions (as confirmed by the agreement between FLIM based inference and OCR measurements in Figures 5 and 8).
Once α is calibrated at the control condition using Equation S48, and βeq is determined from an oxygen drop experiment, then subsequent FLIM measurements of β and [ NADH f ] can be used with Equation S42 and S41 to determine the absolute value of r ox and Jox for all conditions ( Figure 5).
Inferring r ox from NADH long fluorescence lifetime τ l In this section, we derive an alternative procedure for determining the turnover rate of free NADH r ox , and hence Jox , using changes in the NADH long fluorescence lifetime. The NADH long fluorescence lifetime, τ l , is associated with enzyme-bound NADH (Sharick et al., 2018). In the coarse-grained NADH redox model described above, and in Figure 4b, the enzyme-bound NADH consists of reductase-bound NADH ( [ NADH · Re ] ) and oxidase-bound NADH ( [ NADH · Ox ] ). We therefore assume that the experimentally measured NADH long lifetime, τ l , is a linear combination of the lifetimes of [ NADH · Ox ] and [ NADH · Re ] : where τox and τre are the fluorescence lifetimes corresponding to the oxidase-bound NADH and reductase-bound NADH, respectively. Solving for [ NADH · Ox ] and [ NADH · Re ] as a function of β using Equations S39-S40 and S44, and substituting into Equation S49, we predict that the NADH long fluorescence lifetime τ l is linearly related to the inverse of the NADH bound ratio 1/β : The predicted linear relationship between τ l and 1/β is empirically observed during oxygen drop experiments, as shown in Figure 6a and Figure 6-figure supplement 1. This is a self-consistency check that argues for the validity of the assumption in Equation S49.
At equilibrium, when there is no flux through the ETC (i.e., Jox = 0 ), Equation S50 gives: where τeq is the NADH long lifetime at equilibrium. Solving for β and βeq as a function of τ l and τeq from Equations S50 and S53 and substituting into Equation S42, we obtain r ox in terms of τ l : where A and B are the slope and offset of the linear relation between τ l and 1/β in Equation S50. We experimentally measured A and B for each oocyte from the slope and offset of a linear fit between τ l and 1/β during oxygen drop experiments across all drug perturbations (Figure 6a; Figure 6-figure supplement 1). We obtained the equilibrium long lifetime, τeq , by FLIM measurements at the lowest achievable oxygen level in our set up: [ O 2 ] = 0.26 ± 0.04 µM . Once A, B, and τeq are measured, Equation S54 can be used to determine r ox from FLIM measurements of τ l . If α is not known, this procedure can only be used to obtain r ox up to a constant of proportionality. If α is independently measured from Equation S48 at one condition, then Equation S54 can be used to determine the absolute value of r ox for all conditions ( Figure 6b).
As described in the main text, r ox inferred from τ l using Equation S54 produces the same results as r ox inferred from β using Equation S42 (Figure 6b). The agreement between these two methods is a strong self-consistency check of the NADH redox model.

Accounting for NADPH and other background fluorescence
Equation S41 provides a method to infer ETC flux since all factors in it, except for the constant of proportionality, α , depend only on [ NADH b ] and [ NADH f ] , which can be measured from FLIM of NADH. One potential complication with this procedure is that NADPH, another autofluorescent electron carrier, shares a similar fluorescence spectrum with NADH, resulting in a mixed NAD(P) H signal from the autofluorescence measurement. While NADH concentration is 40 times greater than the concentration of NADPH for the whole mouse oocytes (Bustamante et al., 2017), and presumably even higher for mitochondria, NADPH concentration can be comparable to that of NADH for other cell types such as tissue culture cells (Park et al., 2016). In this section, we generalize Equation S41 to predict the ETC flux by explicitly considering the potential contributions of other fluorescence species, such as NADPH, to the measured autofluorescence signal.
We start from the fact that concentrations of bound and free fluorescent species measured from FLIM using Equations S2-S3, [  . If the signal from NADPH and other additional fluorescence species is additive, then: [ where C f and C b are the non-NADH contributions to the measured concentrations of free and bound fluorescent species. We substitute Equations S55-S56 to the predicted ETC flux in Equation S41 and obtain Rearranging, we obtain where Comparing Equation S58 with Equation S41, we notice that the background fluorescence does not change the form of the equation of the predicted ETC flux because the concentrations of the background fluorescent species are incorporated into the equilibrium bound ratio β N,eq . If β N,eq can be reliably measured, the background fluorescence will not affect the flux inference procedures. In other words, can be used for flux inference in place of [ NADH f ] and [ NADH b ] in Equations S41-S42. Therefore, an additive offset to the measured concentrations of free and bound species will not affect the flux inference procedure, whether that additive offset comes from NADPH or from other sources of fluorescent background.
Alternatively, if the signal from background fluorescence changes proportionally with NADH, then: we have where Comparing Equation S63 with Equation S41, we again obtained the same form for Jox , but with a rescaled α N . Therefore, a background fluorescence signal that changes proportionally with NADH will not affect the flux inference procedure, whether that background comes from NADPH or from other sources.
Thus, we have shown that the flux of the generalized NADH redox model reduces to the flux of the reversible Michaelis-Menten model in the limit where the unbinding rate of NADH from the reductase is much faster than any other rates in the model.
We note that the predicted flux-concentration relation for the reversible Michaelis-Menten model (Equation S72) remains exactly the same as the generalized model (Equation S41), but with different expressions for α and βeq as expressed in Equations S69-S70.
Next, we generalize the results in Equation S72 to a detailed model with N oxidase and M reductase, each of which is described by the reversible Michaelis-Menten kinetics. Unpacking the coarse-grained binding and unbinding rates from Equations S24-S25, we obtain where k 1,i and k −1,i denote the binding and unbinding rates of NADH to the ith oxidase.
Connecting the NADH redox model to detailed biophysical models of mitochondrial metabolism In this section, we show that the coarse-grained NADH redox model described above, and in Figure 4b of the main text, can be directly related to detailed biophysical models of mitochondrial