Systems biology identifies preserved integrity but impaired metabolism of mitochondria due to a glycolytic defect in Alzheimer's disease neurons

Abstract Mitochondrial dysfunction is implicated in most neurodegenerative diseases, including Alzheimer's disease (AD). We here combined experimental and computational approaches to investigate mitochondrial health and bioenergetic function in neurons from a double transgenic animal model of AD (PS2APP/B6.152H). Experiments in primary cortical neurons demonstrated that AD neurons had reduced mitochondrial respiratory capacity. Interestingly, the computational model predicted that this mitochondrial bioenergetic phenotype could not be explained by any defect in the mitochondrial respiratory chain (RC), but could be closely resembled by a simulated impairment in the mitochondrial NADH flux. Further computational analysis predicted that such an impairment would reduce levels of mitochondrial NADH, both in the resting state and following pharmacological manipulation of the RC. To validate these predictions, we utilized fluorescence lifetime imaging microscopy (FLIM) and autofluorescence imaging and confirmed that transgenic AD neurons had reduced mitochondrial NAD(P)H levels at rest, and impaired power of mitochondrial NAD(P)H production. Of note, FLIM measurements also highlighted reduced cytosolic NAD(P)H in these cells, and extracellular acidification experiments showed an impaired glycolytic flux. The impaired glycolytic flux was identified to be responsible for the observed mitochondrial hypometabolism, since bypassing glycolysis with pyruvate restored mitochondrial health. This study highlights the benefits of a systems biology approach when investigating complex, nonintuitive molecular processes such as mitochondrial bioenergetics, and indicates that primary cortical neurons from a transgenic AD model have reduced glycolytic flux, leading to reduced cytosolic and mitochondrial NAD(P)H and reduced mitochondrial respiratory capacity.


| Calibration of a flux-based computational model of the mitochondrial respiratory chain
We implemented a previously published (Beard, 2005;Huber et al., 2011) computational model of the mitochondrial RC that incorporates fluxes through the mitochondrial respiratory complexes, ATP production mediated by the F 1 F o ATP synthase, the mitochondrial membrane potential, and nucleotide, ion and proton fluxes across the mitochondrial membranes ( Figure 1a). The model is described in detail in Methods and Supporting Information Appendix S1. We first parameterized the computational model using values from the literature (preferentially from wild-type (WT) primary neurons; see Supporting Information Tables S1-S4 for model description and

| Computational analysis suggests a defect in NADH flux to the respiratory chain
We next utilized the computational model to investigate which molecular alterations could produce the observed respiratory   phenotype. We simulated a standard OCR experiment (sequential addition of oligomycin, FCCP and antimycin A) by partial inhibition of the F1Fo ATP synthase, increase in the H+ leak and partial inhibition of complex III (Supporting Information Appendix S1). We then simulated impairments in individual modules by decreasing the activity of the relevant respiratory complex or increasing the activity of the H + leak flux, and compared the predicted OCR to that obtained by the modelling of "wild-type" conditions (WTsim, no simulated impairments). We considered an OCR metric to be altered if it was predicted to differ from WTsim by more than ±10%, and exceeded an additional threshold calculated using statistical power analysis.
Somewhat surprisingly, the model predicted that defects in complex I (CI), CIII, CIV, F 1 F o ATP synthase (baseline activity reduced to 70% of WTsim) or H + leak (baseline activity increased to 150% of WTsim) would not recapitulate the behaviour observed in TgAD neurons (Figure 2d,e). In contrast, a simulated impairment in the NADH flux to the RC (baseline activity reduced to 95% of WTsim) was predicted to induce a strong decrease in maximal OCR, while leaving the other OCR metrics unchanged (Figure 2d,e,f), thereby recapitulating experimental measurements. Computational analysis therefore suggested that the reduction in maximal respiratory capacity measured in TgAD primary cortical neurons could not be explained by a defect in the RC, but rather by a defect in the provision of NADH (or other substrate) as input to the RC.
We next investigated whether the mitochondrial RC in our TgAD neurons was indeed intact. We measured the protein levels of respiratory complexes I, II, III and IV and the F 1 F o ATP synthase, and indeed did not observe any significant difference between genotypes (Figure 3a,b). As measurements of mitochondrial activity could also be impacted by changes in the numbers of mitochondria or morphology of the network, we labelled mitochondria with a mitochondria-targeted red fluorescent protein (mtDsRed; Figure 3c), and performed quantitative morphological analyses (Koopman et al., 2005). The form factor (FF) and aspect ratio (AR), two geometric indicators of mitochondrial shape in terms of elongation and branching, were similar in WT and TgAD neurons, suggesting unaltered mitochondrial dynamics ( Figure 3d).
Moreover, the size, number and total area of mitochondria were also similar in both genotypes (Figure 3e), collectively demonstrating intact mitochondrial integrity, morphology and network dynamics.
To further assess activity of the respiratory chain, we investigated mitochondrial membrane potential (ΔΨ m ). The computational model predicted that, while impaired respiratory chain activity induces a measurable decrease in baseline ΔΨ m , a defect in NADH supply induces a smaller reduction ( Figure 3f). The model also predicted that such an impairment has a minimal effect on ΔΨ m fold changes in response to oligomycin, rotenone or antimycin A, but has a greater effect on the response to FCCP (Figure 3g). To investigate these predictions experimentally, we measured ΔΨ m using TMRM.
We first measured baseline ΔΨ m by bathing neurons in a saline containing K-gluconate to depolarize the plasma membrane and remove the confounding effect of plasma membrane potential changes on TMRM fluorescence (Tottene, Moretti, & Pietrobon, 1996;Ward et al., 2007). We did not measure a significant difference in baseline F I G U R E 2 Oxygraphy measurements in live cells show impairment of maximal respiration in transgenic AD neurons, and the computational model predicts that this can be explained by a defect in mitochondrial NADH flux. (a) The experimental protocol followed to measure the oxygen consumption rate (OCR) in primary cortical neurons. The classical "mitochondrial stress test" assesses mitochondrial respiratory activity: oligomycin (Oligo, 2 μg/ml) inhibits the F 1 F o ATP synthase, FCCP (0.5 μM) uncouples respiration, and antimycin A (AntiA; 1 μM) inhibits electron flow and abolishes mitochondrial respiration. Nonmitochondrial OCR (OCR remaining after addition of all drugs) is subtracted from all measurements to calculate the displayed OCR metrics (Brand & Nicholls, 2011): basal = mitochondria-specific respiration at rest; ATP synthesis = basal OCR dedicated to ATP production (oligomycin-sensitive respiration); proton (H + ) leak = basal OCR uncoupled from ATP production (oligomycin-insensitive respiration); maximal OCR = OCR upon uncoupling of respiration from ATP production; spare capacity = "spare" OCR available while at rest (maximal-basal). (b) Mean OCR measured in primary cortical neurons from wild-type (WT) and transgenic AD (TgAD) mice. *p = 0.006. n (independent cultures-number of wells): WT 7-37, TgAD 8-43. (c) OCR metrics calculated from measurements in WT and TgAD cortical neurons, as in (a). Coupling efficiency (= ATP synthesis/basal) reports the relative fluxes through the ATP synthase and proton leak pathways. Cell RCR max. and cell RCR basal (= maximal/H + leak and basal/H + leak, respectively) report the efficiency of substrate oxidation in the basal or maximal respiration state. **p < 0.01; ***p < 0.001. n: WT 10-67, TgAD 12-75. (d) Mean OCR metrics predicted by the computational model with no simulated impairment (WTsim), and impairments simulated in the indicated fluxes (100 simulations for each). Impairments were simulated by reducing fluxes to 70% of WTsim for respiratory complex I (CI), CIII, CIV and F 1 F o ATP synthase; increasing the H + leak flux (Hle) to 150% of WTsim; or reducing the NADH flux to 95% WTsim. The dashed lines indicate the statistical thresholds defined to identify whether predicted changes were likely to be measured experimentally (see Section 4). H+ leak and ATP synthesis are reported as oligomycin-insensitive and oligomycin-sensitive respiration, respectively, to allow direct comparison with experiments (see Supporting Information Appendix S1). (e) Heatmap highlighting the OCR metrics that differ between WT and TgAD neurons. The first column illustrates that only maximal respiration differed between WT and TgAD neurons in experimental measurements (shaded dark red). The subsequent columns indicate the changes predicted by the computational model with impairments simulated as in (d). Changes > ±10% compared to WTsim are marked light blue/red, while changes that exceed the statistically defined thresholds are marked dark blue/red. The model predicted that only an impairment in NADH flux could correctly reproduce the experimentally observed behaviour. (f) Predicted OCR metrics in WTsim compared to simulations with impaired NADH flux. The mean value is plotted as an unfilled circle within the boxplots (100 simulations) THEUREY ET AL.

| NAD(P)H autofluorescence measurements validate the computationally predicted mitochondrial NADH defect in transgenic AD neurons
We next utilized the computational model to suggest additional experiments that could further validate the presence of an impaired NADH supply to the RC. The model predicted that such an impairment significantly alters baseline mitochondrial redox status (NADH/ NAD + ) and the redox response to specific pharmacological inhibition ( Figure 4a). While FCCP maximally increases RC activity and hence NADH consumption, this consumption is balanced by an increase in NADH production/import by mitochondria, thereby establishing a steady state. Thus, addition of rotenone (directly inhibiting NADH consumption by complex I) subsequent to FCCP will induce a rapid increase of NADH levels informative of mitochondrial NADH metabolism. In this instance, the model predicted that mitochondrial redox status, following the addition of FCCP plus rotenone, would be markedly reduced by an impairment in the NADH flux ( Figure 4a).
To validate these predictions, we measured NAD(P)H autofluorescence in our primary cultures, using classical epifluorescence microscopy ( Figure 4b). We note that simulated redox status (WTsim) agrees with NAD(P)H autofluorescence measurements in WT cells F I G U R E 3 Primary cortical neurons from transgenic AD mice have preserved respiratory chain protein expression, mitochondrial morphology and network dynamics, while the mitochondrial membrane potential following uncoupling is significantly reduced. (a) Representative Western blots with molecular weight markings to the right of the blots, and (b) corresponding densitometry analysis of the expression level of subunits of the respiratory chain complexes (I-IV) and the F 1 F o ATP synthase, from WT and transgenic AD (TgAD) primary cortical neurons. The heat-shock protein 90 (HSP90) was used as a loading control. n (independent cultures): WT 5 and TgAD 7. Black vertical lines indicate where some blots were cut to remove the molecular weight marker. (c) Representative confocal images of primary cortical neurons (after 6 DIV) transfected with a mitochondrial red fluorescent protein (mtDsRed) highlight the intricate mitochondrial network throughout the neuron. Scale bar, 10 μm.
Total ATP production was also lower in TgAD neurons following oligomycin (−17%, p = 0.003, Figure 6d), further suggesting impaired glycolysis in these neurons. Although total ATP production did not significantly differ between the two genotypes in the basal state ( Figure 6e), the contribution of mitochondrial oxidative phosphorylation to ATP production was nevertheless higher in TgAD neurons, again indicating a more aerobic basal metabolism (+6%, p = 0.001, Figure 6d).
To investigate the source of glycolytic impairment, we measured the expression levels of several rate-limiting enzymes and transporters associated with glucose metabolism-glucose transporter 3 (GLUT3), hexokinase 1 (HK1), lactate dehydrogenase (LDH) and mitochondrial pyruvate carriers 1 and 2 (MPC1-2). Interestingly, we did not observe any differences between the two genotypes ( Figure 6f,g), suggesting that a defect affecting glycolysis might originate from levels of regulation other than the expression levels of these proteins.

| Reduced mitochondrial NAD(P)H in transgenic AD neurons is caused by an impaired glycolytic flux
We further investigated the causal relationship between the glycolytic impairment and mitochondrial NAD(P)H flux reduction. In TgAD neurons, an impaired glycolytic flux, and the resultant reduction in cytosolic NAD(P)H and carbon fluxes to mitochondria, could lead to decreased import/production of NAD(P)H by mitochondria. In this case, we hypothesized that providing pyruvate, rather than glucose, to the neurons would bypass any glycolytic impairment and F I G U R E 6 Extracellular acidification rate (ECAR) measurements demonstrate a glycolytic defect in transgenic AD neurons, and experiments in pyruvate suggest a causal relationship between glycolytic and mitochondrial impairments. (a) Experimental protocol followed to measure the extracellular acidification rate (ECAR) in primary cortical neurons. Maximal glycolysis = ECAR after inhibition of F 1 F o ATP synthase activity with oligomycin (Oligo; 2 μg/ml). Glycolytic reserve (Glyc. reserve) = "spare" ECAR available while at rest (maximal-basal). (b) Basal ECAR and the oxygen consumption rate (OCR)/ECAR ratio in WT and transgenic AD (TgAD) neurons. TgAD neurons had significantly lower ECAR compared to WT, and were more aerobic (higher OCR/ECAR ratio). *p < 0.05, ***p < 0.001. n (independent cultures-number of wells): WT 10-134, TgAD 12-120. (c) Maximal glycolysis (Max. glyc.) and glycolytic reserve were also significantly lower in TgAD transgenic neurons (**p < 0.01, ***p < 0.001). n: WT 5-44, TgAD 7-51. (d) Proportion of glycolytic (black) and oxidative (grey) ATP production in basal conditions. (e) Glycolytic, oxidative and total ATP production in basal conditions and following oligomycin addition (+Oligo; 2 μg/ml), normalized to WT basal levels. (d, e) **p < 0.01; ***p < 0.001. n: WT 10-73, TgAD 12-76. (f) Representative Western blots with molecular weight markings to the right of the blots, and (g) corresponding densitometry analysis of proteins involved in glucose metabolism-glucose transporter 3 (GLUT3; n (independent cultures): WT 5, TgAD 7), hexokinase 1 (HK1), lactate dehydrogenase (LDH), mitochondrial pyruvate carriers 1/2 (MPC1/2) (all n: WT 4, TgAD 4). Actin or HSP90 were used as loading controls. Black vertical lines indicate where some blots were cut to remove the molecular weight marker. (h-k) Measurements in WT and TgAD neurons supplemented with 5 mM pyruvate. (h) Basal OCR, proton leak, ATP synthesis and maximal OCR, as described in Figure 2a  | 11 of 16 restore mitochondrial health. We therefore repeated key experiments in which we observed significant differences compared to WT neurons, supplementing neurons with pyruvate (5 mM) rather than glucose. Interestingly, this ablated all significant NAD(P)H impairments measured in TgAD neurons (Figure 6h-k). As bypassing glycolysis suppressed mitochondrial defects in TgAD neurons, this indicated that a glycolytic defect may be primarily responsible for the mitochondrial NAD(P)H dyshomeostasis observed in these cells.
Both our experimental and computational analyses argued against a direct impairment of the mitochondrial RC in TgAD neurons. Although impairments in RC activity, specifically in complex IV, have been frequently associated with AD (Bosetti et al., 2002;Kish et al., 1992;;Rhein et al., 2009;Yao et al., 2009), findings are not wholly consistent (Choi et al., 2012;Kipanyula et al., 2012), and studies were primarily performed in postmortem AD brains or in situations where clinical pathology was already apparent. In these instances, RC impairments may be related to Aβ toxicity (Pagani & Eckert, 2011). In our TgAD mice (B6.152H), elevated Aβ deposits have been detected in the brain at 3-6 months of age, with decreased complex IV activity and cognitive defects evident at 8 months (Ozmen et al., 2009;;Rhein et al., 2009;Richards et al., 2003). Thus, it is noteworthy that the bioenergetic impairments reported here were identified in primary cortical neurons prepared from pups prior to any increased Aβ production/deposition or the onset of pathological behavioural phenotypes. These results further suggest a mechanism outside amyloid considerations for the bioenergetic impairment described here.
Our work, demonstrating glucose hypometabolism in vitro, echoes clinical research in humans, with glucose metabolism of increasing importance in AD research. Indeed, impaired glucose metabolism is one of the most reliable indicators of disease progression, potentially preceding clinical symptoms by years or even decades (Bubber et al., 2005;Chen & Zhong, 2013). However, the molecular mechanisms linking glycolysis defects to neurodegeneration remain elusive, with contrasting findings on the level of glucose metabolism enzymes and transporters in AD patients (Bigl, Bruckner, Arendt, Bigl, & Eschrich, 1999;Harr, Simonian, & Hyman, 1995), has been observed in AD patients, its induction is protective in cell and animal models of AD (Zhu, Shan, Yuzwa, & Vocadlo, 2014), and inhibition of a glycoside hydrolase prevented cognitive decline, plaque formation and tau aggregates in transgenic mouse models of neurodegeneration (Hastings et al., 2017;Yuzwa et al., 2014). On this aspect, however, studies in appropriate neuronal models are lacking. Moreover, dysregulated Ca 2+ signalling, as previously described in these TgAD neurons (Kipanyula et al., 2012), could also affect glucose/NAD(P)H metabolism (Llorente-Folch et al., 2015;Pancani, Anderson, Porter, & Thibault, 2011), and Wnt signalling was recently suggested as a potential link between neuronal physiology, glucose metabolism and AD (Cisternas & Inestrosa, 2017

| Primary cultures of cortical neurons
Primary neuronal cultures were obtained from cortices dissected from 0-to 1-day newborn mice of either sex, as previously described (Kipanyula et al., 2012). Experiments were performed after 9-11 days in vitro (DIV) unless stated otherwise.

| TMRM and NAD(P)H epifluorescence microscopy
Experiments were performed following the standardized protocols described recently (Connolly et al., 2017). After verifying signal stability, baseline measurements were used to normalize the traces. We utilized a combination of mitochondrial inhibitors to assess the impact on mitochondrial membrane potential (ΔΨ m ) and NAD(P)H, and therefore investigate ΔΨ m and NADH fluxes (Connolly et al., 2017). To measure baseline ΔΨ m , we bathed neurons in a saline containing 130 mM K-gluconate to dissipate the plasma membrane potential (ΔΨ p ) and overcome the confounding effect of baseline ΔΨ p changes on TMRM fluorescence (Tottene et al., 1996;Ward et al., 2007). Drug concentrations are listed in Supporting Information Table S5.

| Oxygraphy, extracellular acidification measurements and ATP calculations
Live OCR and extracellular acidification rate (ECAR) were measured simultaneously using a Seahorse XFe24 Analyzer (Agilent). Saline and equilibration conditions were the same as for live imaging. The classical "mitochondrial stress test" protocol was performed as previously described (Connolly et al., 2017). After calculation of the buffering power of the experimental saline and the proton rate production (Mookerjee & Brand, 2015), mitochondrial and cytosolic (glycolytic) ATP production rates were calculated from OCR and ECAR measurements as described previously (Mookerjee et al., 2017).

| Immunofluorescence and transfection
Transfection of a mitochondria-targeted red fluorescent protein DsRed (mtDsRed), cloned in a pcDNA3 backbone (Filadi et al., 2015), was performed the day after plating, 6 hr before the change from plating medium to growth medium, with Lipofectamine 2000 (Thermo Fisher), following the protocol provided by the manufacturer. Fixation and immunofluorescence were performed as previously described (Filadi et al., 2015). Primary antibodies directed against GFAP and NF200 are detailed in Supporting Information Table S6.

| 2-photon NAD(P)H fluorescence lifetime imaging microscopy (FLIM)
Neurons were equilibrated at 37°C in the same saline as for NAD(P)H autofluorescence experiments, and imaging was performed at 22°C in the same saline with pH adjusted for temperature. FLIM measurements were carried out using the time-correlated single-photon counting (TCSPC) method. Amplitude (A1, A2) and fluorescence lifetime (τ 1 , τ 2 ) parameters were extracted from the decay fitting.

| Protein extraction and Western blotting
For extraction and denaturation of proteins, two different protocols were used: RC complexes, MPC1 and 2 and GLUT3 were extracted in RIPA buffer and denatured at 55°C; HK1 and LDH were extracted in Tris-HCl and denatured at 70°C. After extraction, proteins were quantified using the BCA Protein Assay Kit (Pierce) following the manufacturer's instructions. Detection of proteins was performed by chemiluminescence. The intensity of the bands was analysed using the Fiji software program.

| Morphological characterization of the mitochondrial network
Labelled mitochondria with a mitochondria-targeted red fluorescent protein (mtDsRed) were morphologically analysed as described in Koopman et al. (2005) using the Fiji software (Schindelin et al., 2012).

| Computational modelling
The core computational model was originally developed as a selfcontained thermodynamically balanced model of oxidative phosphorylation and the electron transport chain, based on experimental observations in isolated cardiac mitochondria (Beard, 2005). It has been compared to, and validated by, experimental data from in vitro cardiac and liver isolated mitochondria and in vivo skeletal muscle mitochondria (Bazil, Beard, & Vinnakota, 2016;Dash & Beard, 2008;Wu, Yang, Vinnakota, & Beard, 2007), and was extended (additional flux incorporated to describe cytosolic ATP processes) to analyse respiratory chain function in intact cancer cells (Huber, Connolly, Dussmann, & Prehn, 2012;Huber et al., 2011). We chose to implement this model (Figure 1a) to specifically focus our computational investigations on the mitochondrial respiratory chain, impairments of which have been extensively reported in AD (Bosetti et al., 2002;Kish et al., 1992;Rhein et al., 2009;Yao et al., 2009  Data points outside this range are marked by "+." Differences between medians were determined using a rank sum test (MATLAB, The MathWorks). In all figures, * indicates p < 0.05, ** indicates p < 0.01, and *** indicates p < 0.001. To determine whether differences predicted by simulated impairments might be experimentally measured, we set a threshold using statistical power analysis (dashed lines in Figure 2d; see Supporting Information Appendix S1).