Quantitative modeling of pentose phosphate pathway response to oxidative stress reveals a cooperative regulatory strategy.

Living cells use signaling and regulatory mechanisms to adapt to environmental stresses. Adaptation to oxidative stress involves the regulation of many enzymes in both glycolysis and pentose phosphate pathways (PPP), so as to support PPP-driven NADPH recycling for antioxidant defense. The underlying regulatory logic is investigated by developing a kinetic modeling approach fueled with metabolomics and 13C-fluxomics datasets from human fibroblast cells. Bayesian parameter estimation and phenotypic analysis of models highlight complementary roles for several metabolite-enzyme regulations. Specifically, carbon flux rerouting into PPP involves a tight coordination between the upregulation of G6PD activity concomitant to a decreased NADPH/NADP+ ratio and the differential control of downward and upward glycolytic fluxes through the joint inhibition of PGI and GAPD enzymes. Such functional interplay between distinct regulatory feedbacks promotes efficient detoxification and homeostasis response over a broad range of stress level, but can also explain paradoxical pertubation phenotypes for instance reported for 6PGD modulation in mammalian cells.

Introduction

The oxidative pentose phosphate pathway (oxPPP) is a fundamental pathway of glucose metabolism involved in nucleotide biosynthesis and redox homeostasis (Stincone et al., 2015). Its role is prominent following an oxidative stress to generate NADPH required for fueling antioxidant machinery and producing biosynthetic precursors required for repairing DNA damage. A significant increase of metabolic flux in this pathway is commonly observed in any living cells subjected to oxidative stress (Ben-Yoseph et al., 1996; Ralser et al., 2007; LaMonte et al., 2013; Kuehne et al., 2015; Christodoulou et al., 2018; Nikel et al., 2021). A common rationale for such metabolic flux rerouting relies on the acknowledged roles of NADP as a coenzyme and NADPH as a competitive inhibitor of the first oxidation reaction of the oxPPP (Warburg and Christian, 1936; Negelein and Haas, 1935; Eggleston and Krebs, 1974). The scavenging activity of glutathione antioxidant is coupled to the oxidation of NADPH into NADP+, which is therefore prone to increase G6PD activity and NADPH production. However, oxidative stress has also been shown to induce the allosteric or oxidative inhibition of diverse glycolytic enzymes such as PGI (Kuehne et al., 2015; Dubreuil et al., 2020), GAPD (Ralser et al., 2007; Peralta et al., 2015), PK (Anastasiou et al., 2011), or TPI (Grüning et al., 2011, 2014). A complex pattern of regulation at the levels of both PPP and glycolysis raises the question of their coordination for efficient metabolic rerouting.

The metabolic network combining glycolytic and pentose phosphate pathways displays a complicate branching structure comprising both reversible and irreversible reactions, which obstructs intuitive understanding of multisite metabolic regulation. To investigate complex regulatory patterns, the kinetic modeling framework is often used to disentangle the respective and cooperative roles of multiple feedback regulations (Relógio et al., 2011; Pfeuty et al., 2018; Sander et al., 2019). Several kinetic models have addressed oxPPP dynamics with respect to a specific organism and experimental dataset (Thorburn and Kuchel, 1985; Schuster and Holzhütter, 1995; Kerkhoven et al., 2013). Nowadays, advanced metabolomics and fluxomics studies such as kinetic measurements of concentrations and isotopic labeling patterns provide a rich material to build increasingly reliable and comprehensive kinetic models (Miskovic et al., 2015; Foster et al., 2019; Hameri et al., 2019). Regarding the regulation of the oxidative stress response, such data are available and have already been analyzed in terms of significance and ranking of diverse regulatory hypothesis (Kuehne et al., 2015, 2017; Christodoulou et al., 2018). However, questions remain about the interplay and cooperativity between those different regulatory mechanisms.

In the present study, we aim to build a class of kinetic models inferred from a comprehensive dataset associated with the oxidative stress response of neonatal human skin fibroblasts (Kuehne et al., 2015). On the one hand, the network structure of the model is chosen to fit with the available data and with the objective to understand metabolic flux rerouting for H2O2 detoxification. On the other hand, flux analysis from 13C-labeling data and parameter estimation from flux and concentration data are based on Monte Carlo sampling methods (Saa and Nielsen, 2016; Valderrama-Bahamóndez and Fröhlich, 2019; Theorell and Nöh, 2020) so as generate a representative sample of kinetic models consistent with experimental measurement values and their uncertainties. From such model ensemble, we perform a comprehensive set of analysis regarding parameter distributions, transient dynamical responses, dose-response properties, and gain/loss-of-function phenotypes, which portrays the manner how allosteric and redox regulations contribute to the metabolic response upon oxidative stress. In particular, these analyses converge to the notion that distributed allosteric regulation is required for efficient metabolic rerouting where regulatory mechanisms display both complementary and cooperative roles.

Results

Kinetic modeling of metabolic response to oxidative stress

Kinetic modeling of the early metabolic response to oxidative stress follows the conventional framework for building kinetic models of metabolic pathways (Miskovic et al., 2015; Foster et al., 2019; Hameri et al., 2019). The metabolic and regulatory network structure considered in this study is described in Figures 1A and S1 and the workflow from the integration of metabolomics and fluxomics data to model ensemble analysis is recapitulated in Figure 1B and STAR Methods. The dynamic response of a metabolic network to an oxidative stress perturbation can be described by a set of differential equations derived from the stoichiometry and enzyme-kinetic reaction rates:where represents the concentrations of metabolite species i, denotes the stochiometric matrix, represents the reaction rates associated with enzymes j (see Table S1), is the Heaviside step exposure to extracellular concentrations of H2O2, and represents the enzymatic or regulatory parameters (see Table S2). The steady-state equation is given by where capitalized and denote steady-state concentration and flux vectors.

Figure 1

Modeling workflow

(A) Selected set of metabolic and regulatory pathways involved in oxidative stress response (see also Figure S1). The network comprises the glutathione, glycolytic, and pentose phosphate pathways supplemented with a selected set of allosteric and oxidative regulations (red arrows).

(B) Modeling workflow from the integration of 13C fluxomics and metabolomics data (Kuehne et al., 2015) to parameter estimation of kinetic model ensemble and regulation analysis of transient response, dose response, and gain/loss-of-function phenotypes.

The main specificity of our class of kinetic metabolic model relates to its regulatory scheme. So far, kinetic models of oxidative stress response investigate early intracellular responses at the levels of antioxidant pathways (Adimora et al., 2010; Kembro et al., 2013; Benfeitas et al., 2014) or of glycolytic and pentose phosphate pathways (Thorburn and Kuchel, 1985; Schuster and Holzhütter, 1995; Kerkhoven et al., 2013; Christodoulou et al., 2018), leaving out the TCA cycle. The latter studies either did not include allosteric and oxidative regulations or explored all possible metabolite-enzyme regulations. In contrast, our model specifically addresses the respective roles of a subset of regulatory mechanisms that has been acknowledged to contribute to some respect to rapid metabolic adaptation to oxidative stress, namely NADPH-dependent inhibition of PPP enzymes (Yoshida and Lin, 1973; Holten et al., 1976; Christodoulou et al., 2018), 6PG-dependent inhibition of PGI (Parr, 1956; Kahana et al., 1960; Gaitonde et al., 1989; Kuehne et al., 2015; Dubreuil et al., 2020), oxidative inhibition of GAPD (Ralser et al., 2007, 2009; Peralta et al., 2015; van der Reest et al., 2018), G6P-dependent inhibition of HK, and regulation of other NADPH-consuming or producing reactions (Jeon et al., 2012; Lewis et al., 2014; Fan et al., 2014; Chen et al., 2019). Notable exceptions have been to leave out the inhibitions of PK, TPI, or PFK1 enzymes (Anastasiou et al., 2011; Grüning et al., 2011, 2014; Mullarky and Cantley, 2015) as depicted and explained in Figure S1.

Refined analysis of stress-induced flux redistribution

The stress-induced redistribution pattern of metabolic fluxes can be inferred without any knowledge about kinetic parameters. Although an expected feature of such redistribution is the increased flux in the oxidative branch of PPP, it remains unclear to which extent is such increase and whether the oxPPP flux is rather directed toward nucleotide production or toward the nonoxidative branch of the PPP. To gain a quantitative description of stress-induced redistribution of metabolic fluxes, we reanalyze 13C-labeling data (Kuehne et al., 2015) using a stochastic simulation algorithm-based 13C metabolic flux analysis (SSA-based 13C-MFA described in STAR Methods) to simulate the isotope labeling system (Thommen et al., 2022) and a Monte Carlo sampling to determine posterior distribution of flux parameters. We obtain a distribution of all reaction fluxes (Figures 2A and S2A) associated with an accurate fit of mass isotopomer data (see Figures S2C and S2D). The size of confidence intervals associated with the estimated flux distribution determines whether the flux can be estimated accurately enough to be used as a constraint for kinetic modeling (Figure S2B).

Figure 2

Stress-induced flux redistribution

(A–C) Mean and SD of the distribution of normalized flux rates obtained using stochastic stimulation algorithm (SSA) for 13C-MFA (see STAR Methods details). Estimation is restricted to a set of elementary flux parameters (while other fluxes can be derived from balance equations) where the index i indicates the enzyme and the index j indicates the directionality and where estimation is performed for basal (B, blue) or stress (S, red) conditions. Triangles indicate parameter estimation that is statistically significant based on the relative size of confidence interval distribution (see Figure S2B). (B and C) Representation of the net flux rates in basal (B) and stress (C) conditions.

The estimated flux distribution pattern in absence and presence of oxidative stress can be summarized for the main branches of the metabolic network as depicted in Figures 2B and 2C. The metabolic state in the unstressed condition corresponds to a glycolytic flux mode where a minor fraction () of glucose import flux is diverted toward oxPPP (Figure 2B). Exposure of leads to a significant increase in oxPPP flux up to which is further split between nucleotide production and nonoxidative PPP (Figure 2C). It also leads to a significant reduction of about 3-fold of metabolic flux in the lower glycolytic branch below GAP. This value coincides with the fold-change reduction of PEP concentrations about 2.8 reported in the data (Kuehne et al., 2015). The decrease in PEP concentration, concomitant to that of GAPD flux, justifies the model assumption of not considering oxidative inhibition of PK and allosteric inhibition of TPI mediated by PEP (Figure S1).

In addition to the net fluxes associated with the branching architecture of the metabolic network, it is to note that some directional fluxes could be estimated in the nonoxidative PPP reactions as well as in the reversible PGI reaction, which are valuable information for kinetic model building.

Optimization and inference methods identify a plausible ensemble of kinetic models

Besides the redistribution of metabolic fluxes, oxidative stress response also induces rapid changes in metabolic concentrations below the minute timescale (Kuehne et al., 2015), which together provides a valuable dataset to estimate the parameters of kinetic models described by Equations 1 and 3. Our parameter estimation problem consists in estimating the values of 36 parameters (i.e., 12 parameters are fixed including equilibrium constants to consider thermodynamic constraints (Li et al., 2011)) from a dataset comprising 13 estimated values of fluxes and 12 measured values of concentration ratios at . The procedure combines two classes of global optimization methods, namely an evolutionary genetic algorithm and a Markov Chain Monte Carlo (MCMC) algorithm, following a stepwise strategy recapitulated in Figure 3A and detailed in STAR Methods. First, the (Equation 9 in STAR Methods details) of kinetic models whose parameters are randomly sampled ( runs) lies with a confidence interval between 2.6 and 55.4, thus confirming the need of a parameter optimization procedure to find model parameters consistent with experimental data. Second, an evolutionary genetic algorithm is used as a preliminary step to generate a sample of optimized models whose confidence interval of values lies between 0.8 and 3. Third, a small subset of such local optimum solutions whose is used as initial conditions of MCMC sampling algorithm of the parameter space. The parameter distribution obtained with large enough sampling of accepted steps provides an estimation for parameter uncertainty and defines a statistical ensemble of kinetic model that is called and that is analyzed in details in the following. Such statistical model ensemble is characterized with a confidence interval of between 0.63 and 0.86 for which the estimated distributions of fluxes and concentration ratios fall within the range of experimental uncertainties (Figures 3B and 3C) and whose R-squared values are respectively and . Parameter distributions shown in Figure 3D discriminate between stiff and sloppy parameters for which confidence intervals span from few percents to several order of magnitude of parameter variations. Spectral analysis of correlation matrix confirms indeed the existence of a few poorly estimated parameters generally associated with strong correlation between parameters of a same reaction (see Figures S3A and S3B). Last, a dataset about dynamic and dose-dependent concentration responses has been retained to assess the predictive capability of the plausible set of model (see Figure S3C). The dynamical and dose response are overall well predicted with, respectively, and due to sparse discrepancies such as a lower threshold of 6PG response.

Figure 3

Model selection and parameter estimation

(A) Whisker plots associated with random parameter set (), optimized parameter set using an evolutionary genetic algorithm (), and parameter set () sampled with MCMC algorithm ().

(B and C) Whisker plots of the concentration ratio for the model ensemble (red) as compared with the mean and SD of experimental values (black). (C) Whisker plots of a subset of normalized reaction fluxes for the model ensemble ( in blue and in red) as compared with the mean and SD of estimated values (black).

(D) Violin plots of parameters for the model ensemble where the explored parameter space is represented in white (non-gray).

In summary, the parameter estimation procedure generates a plausible set of kinetic models whose parameters show rather narrow distributions, except some parameters that have a little impact on data adjustment or that can be compensated by the change of other parameters. Importantly, most inhibitory constant parameters (, and ) show a narrow distribution, confirming the role of these regulations in shaping the metabolic response to oxidative stress. In the following, we systematically perform analysis on this model ensemble to draw a statistical picture of the regulatory properties.

Transient metabolic responses following stress display multiphasic time courses

The characteristics of the transient dynamics after oxidative stress exposure and before reaching some steady state give some preliminary insights about the respective contributions of passive and regulated metabolic responses (Figures 4 and S4). The minute time resolution in the time series dataset seemed not sufficient to identify trends arising at the second timescale (Ralser et al., 2009; Christodoulou et al., 2018). In the simulations of the model ensemble , metabolites within a same metabolic module share a similar dynamic response profile (Figure 4). First, PPP metabolites display a rapid and significant monophasic increase. Second, antioxidant NADPH and GSH metabolites display a fast and significant decrease followed by a slower increase. Third, glycolytic metabolites show moderate changes where a fast decrease seems followed by a slower increase. These different temporal response patterns are typically characterized with a biphasic response where a fast passive response to the perturbation is quickly followed by a slower regulated response.

Figure 4

Transient dynamics of stress-induced metabolic responses

Temporal response of model ensemble to a step of H2O2 representing the mean value (blue line), a subsample of 50 trajectories (gray line) and the SD (gray shadow).

(A) Dynamic response of all metabolite species until 5 min.

(B–D) Temporal response of H2O2, G6P, and GAPD activity until 100 s, highlighting multiphasic and rapid adaptation responses.

The biphasic nature of the transient response is illustrated in the case of intracellular concentrations of H2O2 and G6P (Figures 4B and 4C). H2O2 shows a sharp increase about several orders of magnitude whose timescale within seconds relates to the basal degradation timescale . The time profile of H2O2 later displays an overshoot in the time course where decrease of H2O2 mirrors the increase in GSH and NADPH, which coincides with the accumulation of glycolytic metabolites including G6P. G6P shows indeed a rapid decrease until before to increase again to eventually exceed its initial value (Figure 4C). The early decreasing phase of G6P dynamics coincides with an increased consumption through G6PD associated with higher NADP+/NADPH levels while the late increasing phase can only be due to a decreased glycolytic flux through PGI related to the inhibition of glycolysis at the levels of GAPD or PGI. GAPD inactivation occurs indeed very rapidly within few seconds (Figure 4D), being a possible candidate for the accumulation of G6P and detoxification of H2O2 occurring within 10th of seconds.

In summary, the biphasic response observed in simulations of the plausible set of models distinguishes between a detoxification response in less than a second using the reservoir of GSH and NADPH and a metabolic rerouting response from few to tenths of seconds involving the inhibition of glycolysis to quickly restore high G6P levels.

Dose-response analysis identifies rate-limiting reactions

The redistribution pattern of metabolic fluxes is determined for a specific level of H2O2 for which 13C labeling data were available. The oxPPP flux normalized to the glucose import flux, , is around 1 (Figure 2A) which is far below the maximal flux capacity associated with the full inhibitions of GAPD and PRPP enzymes (i.e., ). To explore the detoxification and flux capacity at higher stress level, we perform a dose-response analysis of the kinetic model ensemble focusing respectively on the antioxidant response, the concentration response, and the flux response (Figures 5 and S5). Simulation of the metabolic response at as function of H shows a transition in the detoxification response around (Figure 5A). Below this value, the detoxification activity of GPx increases with oxidative stress level thereby keeping low intracellular levels of H2O2 at the expense of an increased reduced state of glutathione. Above this value, GPx activity saturates such that reaction flux is bounded and cannot increase anymore to compensate for the increase in H2O2 production above some level. Beyond this threshold, H2O2 is eliminated by catalase consistently with the idea of that the rate-limiting enzymes depend on intracellular H2O2 concentrations (Makino et al., 2004; Ng et al., 2007). Another qualitative change of the metabolic response is observed for large enough H2O2 (Figure 5B). 6PG does not quickly reach a steady state and continues to slowly increase, as values differ between 5 min and 30 min. The appearance of a slower equilibration dynamics of 6PG coincides with the saturated kinetics of 6PGD enzymatic reaction associated with increased levels of 6PG below that value of . In parallel to the concentration changes of H2O2 and 6PG, the dose response of metabolic fluxes displays a gradual change from a glycolytic mode to oxPPP mode up to indicating that carbon atoms can be recycled multiple times into oxPPP (Figure 5C).

Figure 5

Dose-dependent profile of stress-induced metabolic response. Response to varying level H of extracellular H2O2 measured at 5 min (full line) and 30 min (dashed line)

(A) Dose response in the glutatione pathway of , GPx, and GSH/GSSG ratio.

(B) Dose response of 6PG metabolite (other metabolites are shown in Figure S5).

(C) Dose response of the fluxes through enzymatic reactions G6PD, TKT, GAPD, and PRPP (normalized to glucose import rate ).

In conclusion, dose-response analysis identifies distinct classes of rate-limiting mechanisms, where saturated activities of antioxidant enzymes restrain H2O2 detoxification for a given NADPH load while saturated activities in PPP enzymes and a leak flux through PRPPs reaction restrain the maximal production rate of NADPH.

Regulation analysis reveals complementary ranges of regulatory efficiencies

The model comprises several regulatory mechanisms likely to contribute to the metabolic flux rerouting to oxidative stress. Some of these regulations directly act at the levels of reactions producing or consuming NADPH notably in oxPPP (, , ) while others operate at the level of glycolysis (, and ). The distribution of estimated values provides a preliminary insight for the need of regulation to reproduce experimental data. Computation of inhibitory strength ( for enzymatic activity divided by two) over the model ensemble reveals that stress condition is associated with significant inhibition of PGI and GAPD () but also significant disinhibition of G6PD (Figure 6A). In contrast, other inhibitory strength is either low or poorly constrained such as the feedback inhibition of HK enzyme.

Figure 6

Regulation analysis based on gain/loss-of-function simulations

(A) Inhibitory strength associated with the inhibition of enzyme j by metabolite i for (blue) and (red). indicates that inhibition reduces enzymatic activity by 2-fold.

(B) Sensitivity factor (Equation 10) of output Y with respect to several deletions of regulation in abscisse, where the outputs Y are the metabolic responses of H2O2, NADPH, and G6P to oxidative stress relative to basal levels. Bars are mean values and error bars are standard deviations over the kinetic model ensemble .

(C) Sensitivity factor of H2O2 with respect to deleted regulation , , as function of H. The dose-specific areas of regulatory effect associated with each deleted regulation () is shown upper to the panel.

(D) Response coefficients and (Equation 11) associated with modulated inhibition of upper and lower glycolysis and as function of H.

(E) Second-order sensitivity factor with respect to the combined deletion of two regulations in comparison to the single deletions.

A more comprehensive strategy to quantify regulatory effects consists in measuring gain of function or loss of function associated with the deletion of a single regulatory mechanism (), other things being equal. For such aim, we define a sensitivity quantity measuring the impact of modulating regulatory parameters on some functional output such as intracellular H2O2 (Equation 10 in STAR Methods). larger or lower to one indicates that deleting the regulation increases or decreases the steady-state level of Y. In Figure 6B, removal of NADPH-dependent inhibition of G6PD and H2O2-dependent inhibition of GAPD leads to higher H2O2 and lower NADPH, highlighting significant contributions of these regulations for NADPH homeostasis and H2O2 detoxification. Surprisingly, deleting the inhibition of PGI has a minor impact on metabolic outputs of interest, which can nevertheless be interpreted by a low net flux through PGI () for a particular stress level (Figure 2C).

To investigate whether regulation efficiency indeed depends on oxidative stress level, we evaluate the deleterious effect of removing a regulation (i.e., loss of function) as function of H (Figures 6C, 6D, and S6). The dose-dependent profile of distinguishes different ranges of stress level for which each regulation is the most efficient (Figure 6C). First, the efficiency of G6PD upregulation is the highest for low-to-moderate stress level because enzyme activity rather than G6P level is rate-limiting. Second, the efficiency of 6PG-dependent inhibition of PGI (null for in Figure 6B) peaks at intermediate stress level as it requires both high enough stress for 6PG accumulation (see Figure 5B) and not-too-high stress such that the downward flux prevails over the upward flux in PGI reaction ( or equivalently ) (see Figure 5C). Third, the efficiency of GAPD inhibition culminates at high oxidative stress consistently with its ability, in addition to restore G6P levels, to reverse glycolytic flux coming from the nonoxidative branch of the PPP, thus making possible a flux cycling mode where .

To clarify how to reconcile the results that inactivation of GAPD starts at low intracellular oxidant levels (i.e., ) and GAPD inactivation is the most effective for high oxidative stress (), we also compute response coefficients (Equation 11 in STAR Methods) as function of stress level H (Figure 6D). The results show that the effect of on is already effective for low H values but starts to increase for higher stress level. This dose-dependent effect might relate to the fact that GAPD inhibition is more prone to promote a glycolytic flux toward G6P and oxPPP flux when the carbon flux in the nonoxidative branch of PPP is already flowing from R5P to GAP (due to increased levels of R(u)5P metabolites).

Besides complementary efficiency ranges, two regulations can also combine their effect in non-trivial manners. For instance, G6PD upregulation leads to 6PG accumulation which potentializes the inhibition of PGI. Alternatively, PGI inhibition reduces G6P consumption which potentializes the upregulation of G6PD. These are second-order effects that are quantified by the computation of sensitivity factors associated with the combined deletion of two regulations (Figure 6E).

Ambivalent role of 6PGD enzyme in oxidative stress response

G6PD, 6PGD, and TKT enzymes are common targets of loss-of-function experiments to investigate oxidative stress response (Kuehne et al., 2015; Nóbrega-Pereira et al., 2016; Wan et al., 2017; Sun et al., 2019; Dubreuil et al., 2020). We therefore perform simulations of our model ensemble while varying the enzymatic activity parameter from 10-fold reduction to 10-fold increase. We record both the mean of intracellular H2O2, the NADPH-producing oxPPP flux, and some concentration metabolites for the model ensemble (Figure 7). Expectedly, increase (resp., decrease) of leads to a more (resp., less) efficient oxidative stress response related to subsequent change in the oxPPP flux of NADPH production (Figures 7A and 7B). In sharp contrast, modulation of leads to a more surprising and ambivalent metabolic response (Figures 7C and 7D). Depending on the level of stress and of modulation, we observe that both the increase and the decrease of can hamper the oxidative stress response, while decreasing can both weaken and improve oxidative stress response. These ambivalent phenotypes relate to the dual effect of modulating on G6PD activity itself but also on the increased level of 6PG that inhibits PGI activity and upregulates both G6P level and G6PD flux. Specifically, reducing can lead to an imbalanced state associated with increased metabolic flux through 6PGD and decreased flux through G6PD, whose relative effect determines the occurrence of loss or gain of function. Last, the modulation of has a moderate effect on flux reprogramming restricted to higher stress level (Figures 7E and 7F), consistently with the dose-dependent increase of flux in nonoxidative PPP (Figure 5C) and its crucial role for enabling a cycling flux where . It is to note that some trends observed for the concentration response of glycolytic and PPP intermediates in response to modulation of and (see Figure S7) are qualitatively consistent with knockdown experiments (Kuehne et al., 2015).

Figure 7

Gain/loss of function associated with modulated activity of PPP enzymes

Sensitivity of output Y with respect to modulation of X (Equation 10) as a function of H and the extent of parameter modulation.

(A–F) The output variable Y is either H2O2 (A, C, and E) or NADPH-producing flux (B, D, and F) and the modulated parameter is either (A and B), (C and D) or (E and F).

Discussion

In this work, we use a kinetic modeling approach to investigate the interplay between several metabolic regulations in that specific functional context of oxidative stress response. Data-driven kinetic modeling framework is commonly used to identify large-scale patterns of metabolite-enzyme interactions and to investigate the role of allosteric regulation in controlling metabolic phenotypes (Grimbs et al., 2007; Link et al., 2013; Machado et al., 2015; Jahan et al., 2016; Reznik et al., 2017; Millard et al., 2017; Christodoulou et al., 2018; Britton et al., 2020). The present study purposely focuses on a selected set of regulation proposed to contribute to the rapid and global changes of flux and concentration upon oxidative stress, resulting in PPP upregulation and NADPH recycling to support antioxidant systems (Ralser et al., 2007; Kuehne et al., 2015; Mullarky and Cantley, 2015; Stincone et al., 2015; Dick and Ralser, 2015; Christodoulou et al., 2018). By integrating metabolomics and fluxomics data, kinetic modeling reveals that the observed concentration and flux responses require the co-operation of multiple regulatory mechanisms, and that such co-operation follows a rationale of metabolic efficiency in terms of complementary and synergistic effects. Specifically, we highlight that the regulations of several PPP and glycolytic enzymes all fulfill different but coordinated roles for stress-induced carbon flux rerouting.

Using NADP+ as a substrate and NADPH as a feedback inhibitor, G6PD enzyme is sensitively activated by the decrease of NADP+/NADPH redox ratio associated with the oxidation of NADPH (Yoshida and Lin, 1973; Holten et al., 1976; Christodoulou et al., 2018). However, the efficiency of this sole regulatory mechanism is shown to be hampered by the concomitant decrease of G6P when oxPPP fluxes becomes of the order of downward glycolytic flux. Mechanistically, such decrease of G6P can be compensated either by a reduction of G6P consumption through PGI inhibition or an increase of G6P production through accumulation of glycolytic intermediates itself mediated by the inhibition of glycolytic enzymes below F6P. In our class of kinetic models built from a specific dataset, such inhibition is likely to be mediated by the oxidative inactivation of GAPD, but the inhibition of TPI (Grüning et al., 2011) or PFK1 (Seo and Lee, 2014) could also fulfill a similar flux rerouting role. Interestingly, inhibition of GAPD is more effective to promote oxPPP for significant stress levels coinciding with a carbon flux in the nonoxidative branch of PPP flowing from accumulated R5P to glycolytic intermediates. The inhibitory effect of a given glycolytic enzyme tightly depends on its coupling with nonoxidative PPP reactions at the level of F6P and GAP metabolites. Each regulation mechanism thus implements a specific flux rerouting strategy which results in complementary ranges of efficiency but also synergistic cooperative effects, providing altogether a flexible metabolic adaptation for diverse stress and cellular contexts. A systematic mapping between the regulation motifs and the flux control properties (Machado et al., 2015; Britton et al., 2020) is needed to refine our understanding of the role of regulation for maximizing NADPH yield over a broad range of perturbation and irrespective to the metabolic state.

Besides insights into the regulatory logic, the modeling approach could also identify the main rate-limiting processes for PPP flux increase and H2O2 detoxification. In the models, the maximal PPP flux generates up to 3 mol of NADPH for 1 mol of glucose due to a cycling mode where each molecule of glucose can be oxidized multiple times (Kuehne et al., 2015; Dick and Ralser, 2015; Britt et al., 2022), which is nevertheless far from the maximal yield of 12 mol of NADPH in case of complete inhibition of PRPP and GAPD enzymes. Besides the production of nucleotide precursors for DNA damage repair and of ATP for stress management, we also identify several rate-limiting reactions including GPx enzyme (Ng et al., 2007), but also 6PGD enzyme. For instance, the tight regulation of 6PGD activity is critical for a functional accumulation of G6P. Indeed, both upregulation and downregulation of 6PGD can penalize the oxidative stress response by reducing PGI activity, either too weakly for enabling glycolysis shunt at low stress and too strongly for enabling carbon recycling at high stress. Paradoxically also, the inhibition of 6PDG can both lower the maximal flux capacity for NADPH production and promote the shunt of glucose into oxPPP through the inhibition of PGI. This ambivalent effect explains the contradicting experimental observations where genetic or pharmacological inhibition of 6PGD can either promote or decrease oxidative stress response of mammalian cells depending on the context (Sun et al., 2019; Liu et al., 2019; Dubreuil et al., 2020). Knowing that AMPK, a master regulator of many metabolic processes, is itself regulated by metabolites in the oxPPP (Lin et al., 2015; Gao et al., 2019) and regulates glycolysis especially PFK1 (Wu and Wei, 2012), a better understanding of the regulated coordination between glycolysis and PPP is important for a broad range of metabolic adaptations beyond oxidative stress response.

Limitations of the study

In the present study, kinetic models have been built from a specific experimental dataset associated with neonatal human skin fibroblasts exposed to H2O2-supplemented medium, which questions whether unveiled regulatory patterns can be extrapolated to other types of cells and other sources of oxidative stress. Although similar metabolic responses have been observed for different mammalian cell types and oxidative stress sources (Kuehne et al., 2015; van der Reest et al., 2018; Christodoulou et al., 2019), there are also evidences of alternative metabolic adaptation strategies which restricts the generalization of our findings. For instance, the feedback loop involving the regulation of PK and TPI glycolytic enzymes has also been shown to contribute in oxidative stress response, especially in yeast (Anastasiou et al., 2011; Grüning et al., 2011, 2014). Alternatively, oxidative stress response in some lower eukaryotic and prokaryotic organisms has been shown to involve the Entner-Doudoroff pathway (Nikel et al., 2021) or the ribose salvage pathway (Xu et al., 2013) which are both directly coupled to the PPP. Another limitation of the study inherent to any kinetic model confronted to parameter identifiability issues relates to the possible biases related to model abstractions and assumptions about the choice of reaction kinetics and of involved metabolic pathways leaving aside for instance large-scale metabolic activities and adaptations in the metabolism of ATP, NAD(H), NADP(H), or glutathione.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the Lead Contact, Benjamin Pfeuty (benjamin.pfeuty@univ-lille.fr).

Materials availability

This study did not generate new unique reagents.

Method details

SSA-based 13C-MFA

Stochastic simulation algorithm (SSA) for 13C-based metabolic flux analysis (13C-MFA) is a direct method for the forward simulation problem to compute the dynamics and steady state of (mass) isotopomer distribution in isotopic labeling networks (Thommen et al., 2022). From a given flux distribution and initial labeling state, SSA computes the temporal evolution of isotopomer numbers which are pooled to obtain the mass isotopomer distribution of species j. The sample size parameter of the algorithm is . From the experimental values measured for mass isotopomer, species j in labeling conditions (Kuehne et al., 2015), SSA is performed iteratively using a MCMC sampling method based on a random walk Metropolis algorithm (see section quantification and statistical analysis) for obtaining the posterior distribution of flux parameters. In such flux estimation procedure, the root mean square error function is given by:where the measurement time is and the number of experiments is (, , ).

Differential equation model

The temporal behavior of the metabolic network is described by a differential equation system:

For reversible reactions, we also introduce directional fluxes satisfying the relation where, by convention, the direction goes from G6P to GAP in glycolysis and from R5P to GAP in nonoxidative PPP. Note that represents the net flux of the two reactions between F6P and FBP, where corresponds to the PFK1 reaction and corresponds to the FBPase reaction.

List of abbreviations for metabolites and enzymes

Glc, intracellular glucose; G6P, glucose-6-phosphate; F6P, fructose-6-phosphate; FBP, fructose-1,6-bisphosphate; ALD, fructose 1,6 bisphosphate aldolase; DHAP, dihydroxyacetone phosphate; GAP, glyceraldehyde-3-phosphate; PEP, phosphoenolpyruvate; 6PG, 6-phosphogluconate; 6PGL, 6-phosphogluconolactone; X5P, xylulose 5-phosphate; R5P, ribose 5-phosphate; TPI, triosephosphate isomerase; Ru5P, ribulose 5-phosphate; OX, oxidative stress E4P, erythrose-4-phosphate; S7P, sedoheptulose 7-phosphate; PEP, phosphoenolpyruvate; GT, total glutathione; NT, total NADP; Cat, catalase GR, glutathione reductase; GPx, glutathione peroxydase; Prx, peroxiredoxin; Trx, thioredoxin; HK, hexokinase; G6PD, G6P dehydrogenase; 6PGD, 6PG dehydrogenase; GLase, 6-phosphogluconolactonase; PRPP, phosphoribosyl pyrophosphate; PGI, phosphoglucose isomerase; PFK, phosphofructokinase (type 1); FBPase, fructose-1,6-bisphosphatase; GAPD, GAP dehydrogenase; PK, pyruvate kinase; TLD, transaldolase; TKT(1/2), transketolase (type 1/2); NNH, conversion from NADP+ to NADPH; NHN, conversion from NADPH to NADP+.

Enzyme-kinetic reaction rates

Reaction rates can be described by diverse kinetic laws (e.g., mass action, Michaelis-Menten, Hill, Monod-Wyman-Changeux…). In the present work, we mainly use the generalized mass action kinetics, eventually associated with inhibition:where m and n are the number of substrates and products.

However in the context of oxidative stress response, the significant increase of some metabolite concentrations may require to take into account saturation constants that appear due to the formation of intermediate complexes. A generalization of the uni-uni Michaelis–Menten equation can also be used for irreversible and reversible uni/bi-substrate reaction kinetics (Rohwer et al., 2006):where for unisubstrate reactions and for bisubstrate reactions, and where the denominator contains numerous terms for each substrates and products. In some bisubstrate reactions, the assumption that is made so as to reducing the number of terms and parameter in denominator. Specifically, Michaelis constants are considered only for 6PG, H2O2, GSSG and GSH and in nonoxidative PPP reactions. For GPx and GR reactions, the saturation terms follow a common scheme based on experimental data (Benfeitas et al., 2014). A minimal description of saturation in nonoxidative PPP reactions consists in keeping only product terms (which prevails in case where larger than one) and assumes a same Michaelis constant for substrates and products.

The oxidative stress-dependent inhibition of NADPH consumption related to putative biosynthetic shutdown (Fan et al., 2014) is modeled as an effective rate constant . As well, a simplified description of the oxidative inhibition of GADP enzyme considers a two-state enzyme (non-oxidized and oxidized) where the back-and-forth transition rates depend on the concentration of H2O2 and GSH, respectively:which becomes following quasi-steady-state approximation and linear approximation of regulatory functions and :from which we derive the expression of given in Table S1. The reaction kinetic law associated to each enzymatic reaction is recapitulated in Table S1.

Parameter setting and estimation

Model includes 48 parameters including kinetic constants , equilibrium constants , Michaelis (saturation) constants , inhibitory constants , where i denotes the enzyme, but also conserved quantities and . To restrict and ease the parameter estimation process, parameters are bounded within a restricted range of realistic values as it is often assumed (Christodoulou et al., 2018). An exception is that 8 equilibrium constant values have been fixed based on estimated physiological standard free energy of corresponding reactions (Li et al., 2011). In addition, two saturation constants and the conserved quantities and have also been fixed based on typical values (Benfeitas et al., 2014) and the modeling requirement to reduce parameter space and to constraint parameter estimation. A specific model of index j is denoted as and a set k of model is denoted .

Given a parameter space where parameter values are restricted by upper and lower bounds (listed in Table S2), the parameter estimation problem consists in scoring parameter set by using a normalized root mean square error (nRMSE) :where are concentration ratios in log scale and are the estimated fluxes of reaction i in basal condition and in stress condition. Experimental data consists in and values of concentration ratio and estimated fluxes with a typical standard deviation and that is estimated from 13C-MFA.

This error function is used first with a population-based metaheuristic algorithm called evolution strategy to quickly explore and sample local solutions (Pfeuty and Thommen, 2016). Starting with a pool of parameter vectors of random values (bounded uniform distribution), the algorithm involves three steps: (i) a reproduction step where a parent is randomly selected to be duplicated without applying any fitness criterion at this stage, and to generate offspring; (ii) a mutation step where the parameters of each offspring are modified with a probability through multiplication by a factor , where is a random number of uniform distribution; (iii) a selection step where the nRMSE of the offspring are evaluated, and only the highest-fitness individuals in the pool of parameter sets (parents and offspring) are selected to generate the parents of the next generation. Finally, the optimization process terminates after a maximal number of generations () where a local minimum of is usually reached. This evolution is repeated for a broad set of random initial conditions to obtain a first set of optimized models, some of which satisfies the plausibility criteria . The lowest-nRMSE model is used as an initial condition of the MCMC algorithm (see section quantification and statistical analysis) used to sample a posterior probabilistic distribution of parameters. After a transient, a sample of accepted values are recorded to converge to a representative model ensemble associated to quasi-stationary distributions, while a random subsample of models is used for the statistical analysis of model behaviors.

Sensitivity analysis

Sensitivity quantities are defined to evaluate the effect of parameter changes on stress response behaviors. A model ensemble is defined by a set of model . For a given kinetic model of parameter values , a parameter is multiplied by a factor . A sensitivity factor that informs about the impact of a given parameter on the output state variable Y in response to a stress level H is given by:where . The choice of parameter and variation are selected (i) to evaluate the impact of a regulatory mechanism (i.e., are inhibitory parameters and ) and (ii) to mimmic overexpression or knockdown experiments (i.e., are enzymatic activity rates and ). Another class of sensitivity quantities are response coefficients:

The relation between both sensitivity quantities is given by .

Model updating procedure

From a methodological perspective, the modeling approach relies on a systematic procedure for model inference and parameter estimation, which deliberately considers more parameters than data to provide the possibility to update the probability distribution as additional data becomes available or as additional mechanistic hypothesis are tested. MCMC method is very well adapted for such model updating procedure (Beck and Au, 2002). Specifically, MCMC sampling can be restarted from a given parameter distribution with a modification of the score function and of the sample space. By testing the predictive capacity of the model on dataset that are not included in the model inference procedure, we could indeed identify sparse phenotypes poorly described by the model. An example of model improvement would be to add the corresponding dataset into the monte-carlo sampling procedure which might slightly shift the parameter set of solution to be consistent with these new dataset.

Quantification and statistical analysis

Monte Carlo Markov Chain (MCMC) method

Monte Carlo Markov Chain (MCMC) method is used for estimating the distribution of (i) metabolic flux from 13C labeling data and (ii) kinetic model parameters from concentration data and flux estimation. These two computational problems are both defined by a set of parameters , a set of data values and a root mean square error measuring the difference between data predicted by the parameters and the actual data (see Equations 2 and 9). MCMC simulation methods use Monte Carlo sampling techniques to build Markov chains that converge to the posterior distribution of parameters θ associated to data Y. For both problems, we consider a bounded uniform distribution for the prior. For sampling, a Metropolis-Hasting algorithm is used to generate the random walk Markov chain based on (i) a Gaussian jumping distribution and (ii) an acceptance rate function α that is the ratio of likelihood associated to the next parameter state and actual parameter state. Typically, the likelihood function can be written based on the expected or assumed distribution of observed data, as . Specifically, we use the acceptance rate where the variances of jumping distribution and likelihood function σ are chosen to provide reasonable acceptance rate () and convergence rate. a.

R-squared

Although R-squared is commonly used to assess the goodness-of-fit for a regression model, this quantity can be applied, in addition to RMSE, to assess the fitting and prediction score for the ordinary differential equation models without requiring knowledge about experimental uncertainties:where and Y represents some measurement such as fold-change concentrations.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
Target experimental data	https://doi.org/10.1016/j.molcel.2015.06.017	Figure 2; Tables S1 and S5
Matlab code and simulation data	https://doi.org/10.5281/zenodo.6656078https://github.com/JHurb/HurbainPaper_Algorithm/tree/v1.0.0	JHurb/HurbainPaper_Algorithm-v1.0.0.zip
SBML model	http://www.ebi.ac.uk/biomodels/MODEL2204220001	MODEL2204220001
Software and algorithms
Matlab	https://www.mathworks.com/products/matlab.html	R2018a
SSA 13CMFA	arXiv:2201.00663v1