Control Analysis of Cooperativity and Complementarity in Metabolic Regulations: The Case of NADPH Homeostasis

Complex feedback regulation patterns shape the cellular metabolic response to external or internal perturbations. We propose here a framework consisting of a sampling-based metabolic control analysis of kinetic models to investigate the modes of regulatory interplay in metabolic functions. NADPH homeostasis, for instance in a context of oxidative stress, is an example of metabolic function that involves multiple feedback regulations which raises the issue of their concerted action. Our computational framework allows us to characterize both respective and combined effects of regulations, distinguishing between synergistic versus complementary modes of regulatory crosstalk. Synergistic regulation of G6PD enzymes and PGI enzymes is mediated by congruent effects between concentration sensitivities and reaction elasticities. Complementary regulation of pentose phosphate pathway and lower glycolysis relates to metabolic state-dependent range of regulation efficiency. These cooperative effects are shown to significantly improve metabolic flux response to support NADPH homeostasis, providing a rationale for the complex feedback regulation pattern at work.


Introduction
Metabolic control analysis (MCA) provides a rigorous theoretical framework to study the sensitivity of metabolic networks with respect to biochemical and environmental variations [1,2]. MCA has been further developed and expanded in several directions related to regulation, thermodynamics, or statistical analysis [3][4][5][6][7][8][9]. These developments contribute to more comprehensive analysis of control properties of metabolic networks with the challenging goal to decipher the logic of complex regulation pattern, such as those involving direct metabolite-enzyme interactions and coupling distal parts of a network [10][11][12].
The role of complex feedback regulatory scheme in shaping the network response to environmental changes is recognized in many contexts ranging from nutrient utilization, end-products homeostasis, or stress response [13][14][15]. This is illustrated by the cellular function of NADPH homeostasis, which involves the concerted action of a broad set of metabolic regulation [16][17][18]. NADP(H) homeostasis, like NAD(H) homeostasis [19], is important to keep a functional redox balance against diverse perturbations due for instance to oxidative stress, metabolic stress [20], or reductive biosynthesis [21]. Among the few metabolic pathways producing NADPH [22,23], the oxidative branch of the pentose phosphate pathway (oxPPP) is the main source of NADPH and is also stringently regulated by a set of allosteric and oxidative regulations [16]. Since the characterization of the feedback inhibition of G6PD and 6PGD through competitive binding of NADPH [24][25][26], the metabolic regulatory picture has became increasingly refined with evidences of a significant role for 6PG-dependent inhibition of PGI [17,27] or for oxidative inhibition of several glycolytic enzymes such as GAPD, PFKFB3, or PKM2 [28][29][30][31]. The regulation of NADPH homeostasis is therefore a valuable case study to assess methods and framework to investigate complex regulation patterns.
To address systemic properties of regulation patterns, our kinetic modeling framework extrapolates metabolic control analysis beyond a reference state, by combining mathematical analysis of control equations and sampling analysis of kinetic space. The idea is to represent the statistical distribution of control coefficients on low-dimensional subspaces defined and constrained by control equations. We first depict a global picture of the control pattern related to NADPH homeostasis driven by oxPPP in the absence of regulation, which reveals some trends which are in contradiction with experimental evidences. Combined mathematical and sampling analysis of control pattern further reveals how the presence of feedback regulation promotes PPP flux rerouting and NADPH homeostasis, involving also synergistic and complementary modes of cooperation. NADPH-dependent inhibition of G6PD and 6PG-dependent inhibition of PGI exhibits synergistic cooperation due to congruent concentration control of 6PG and G6P metabolites. Such a regulatory scheme must be supplemented with feedback inhibition of lower glycolysis to extend the efficiency range of NADPH homeostasis beyond a reference flux state.

Materials and Methods
To investigate metabolic control properties associated with feedback regulations, we consider the following kinetic modeling framework restricted to the main pathways involved in PPP-driven maintenance of NADPH homeostasis ( Figure 1A). The concentration and flux dynamics of metabolic networks is commonly described by an ordinary differential equation system (see Appendix A), where vectors and matrices are denoted with bold and italic font styles, respectively. N defines the stochiometric matrix, s i=1,n m the concentration of metabolite species, and v i=1,n r the reaction rate functions described by mass-action kinetics involving rate constants k i and equilibrium constants K i gathered in p = {k i , K j } while regulatory parameters r i are treated separately. Steady-state concentrations S satisfy: N v(S(p, r), p, r) = 0, and the associated steady-state flux vector is noted J(p) = v(S(p), p), where capitalized letters indicate steady-state quantities. Implicit differentiation of the steady-state Equation (2) with respect to kinetic parameters establishes a matrix expression for control coefficients and elasticities [32]: where the control and elasticity matrix coefficients are given by C (2) and (3), we developed a strategy to investigate the role of regulation in determining the range of variation (lower/upper bounds) and state-dependency of the control coefficient C * of interest (here C * = C S,nh i for NADPH homeostasis and

Given Equations
The general form of the function F and the asymptotic behaviors of such function for some small or large values of r i , S i or J i provides key information about the manner how regulation impacts the control coefficient of interest (e.g., promotes NADPH homeostasis). These functions define low-dimensional manifolds in the space of control coefficients and steady-state variables. The distribution of control coefficients obtained from random sampling of model parameters p can be represented on such manifold (C * (p) = F (C S,J (p), S(p), J(p), r)) to better visualize and discriminate the regulatory and context-dependent features determining control properties.  (4)) and sampling analysis of regulatory crosstalk.

Distribution of Control Coefficients in Absence of Feedback Regulation
To provide a preliminary statistical picture of control patterns associated with NADPH homeostasis, a distribution of control coefficients can be obtained from random sampling of kinetic parameters in the absence of feedback regulations ( Figure 2). This sampling approach allows us to distinguish between different contributions in the variability of control pattern: the variability due to changes in the flux state J or due to the difference in parameters associated with a given flux state J. Parameter sampling is therefore subdivided into k parameter subsets p k (size > 10 4 ) associated with specific flux states J k in the space of elementary flux modes ( Figure 2A). The sampled flux state is restricted to the glucoseconsuming modes of the PPP, which corresponds for the kinetic model of Equation (A1) to a two-dimensional triangle polytope as function of J gapd /J hk and J g6pd /6/J hk . In this space, the glycolytic, nucleotide-producing, NADPH-producing modes are the extremity of such polytope, and we compare control distribution in well distinct domains inside the polytope.
Parameter sampling confirms expected trends in the control pattern with many signdefinite coefficients ( Figure 2B,C). NADP + binding to G6PD provides a primary source of oxPPP flux increase in response to NADPH depletion where 0 < C J,ppp gr < 1 (and −1 < C S,nh gr < 0) without the need of regulations. In addition, the main flux-controlling steps are G6PD, but also PFK1 and PGI, which is consistent with the notion that reduced enzyme activity in upper glycolysis leads to flux rerouting into oxPPP. However, some other features of the control pattern do not seem to match with expectations or experimental evidences. For instance, C J,ppp gr and C S,nh gr significantly decrease for flux state domains characterized with J pgi < 0, indicating an unlikely context-dependency of NADPH homeostasis.
As well, 6PG and G6P concentration control coefficients C S,6pg/g6p gr are strongly negative contradicting the numerous experimental evidences reporting a few-fold increase in 6PG and a moderate increase in G6P in response to oxidative stress [17,27,31]. These features suggest the involvement of regulation to enable C S,6pg/g6p gr > 0 and to increase C J,ppp or C S,nh for a broader range of flux states.

Feedback Inhibitions of PPP and Upper Glycolysis Synergistically Cooperate for Efficient PPP Flux Rerouting
From the matrix expression relating control and elasticity coefficients, one can derive a control equation for NADPH homeostasis as function of a set of regulatory parameters. Keeping in mind that oxPPP flux control and NADPH concentration control are related by C S,nh gr = C J,ppp gr − 1 (Equation (A8a)), we can derive in the Appendix C.1 the following general equation: where elasticities i (see Equation (A5)) depends on regulatory parameters r i . This equation is very informative and can be analyzed in several asymptotic limits (Equation (A11)) to dissect how C J,ppp gr depends on regulatory crosstalk, while a sampling approach exploring the parameter space of the kinetic model is required to confirm or refine the results beyond particular assumptions ( Figure 3).
In the absence of regulation, an increased consumption of NADPH enhances the PPP flux control due to the concomitant increase in NADP + as a cofactor of G6PD enzyme, which provides a maximum flux control of max(C J,ppp gr ) = S nh (maximum for J ppp = 0) (Equation (A11a) and Figure 3A). The regulation r 1 (NADPH-dependent inhibition of G6PD) can efficiently promote such PPP flux control to a maximum extent of C J,ppp gr /S nh = 1 + r 1 for small enough S nh and J ppp (Equation (A11b) and Figure 3B). The regulation r 3 alone (6PG-dependent inhibition of PGI) does not promote PPP flux control (Equation (A11d) and Figure 3C). In sharp contrast, such allosteric regulation strongly enhances PPP flux control in presence of r 1 (Equations (A11e) and (A11f) and Figure 3D). This synergistic effect coincides with a positive control of G6P, which itself requires a strong positive control of 6PG mediated by r 1 (low panels of Figure 3B,D). The importance of a positive concentration control of 6PG and G6P is confirmed by the loss of synergistic effect for r 2 ≥ r 1 related to a loss of positive concentration control for G6P and 6PG ( Figure 3E).  (5) and (A8) (colored arrows) recapitulating the interplay of r 1 , r 2 , and r 3 and key steady-state variables on the control associated with NADPH homeostasis.
In the general control Equation (5), C J,ppp gr is expected to decrease with J ppp while the effect of 2 requires high glycolytic flux J pgi > 0. We therefore check that the synergistic interplay between r 1 and r 3 in promoting PPP flux control is indeed compromised for increasing (resp., decreasing) values of J ppp (resp., J pgi ) ( Figure 3F). Finally, the key roles of S nh and C 6pg can be depicted by plotting the distribution of PPP flux control of sampled model on a surface derived from Equation (A9) ( Figure 3G). To summarize, the role of regulations in shaping NADPH homeostasis can be schematically represented to make apparent regulatory crosstalk and context dependencies ( Figure 3H).

Ros-Dependent Inhibition of Glycolytic Enzymes Expands NADPH Homeostatic Abilities
After identifying a synergistic mode of allosteric regulation which is only efficient when glycolytic flux out competes oxPPP flux (J pgi > J ppp ), we now examine the requirement for alternative regulatory strategies in the case where J ppp > J pgi , such as during acute oxidative stress. Indeed, an excessive, endogenous, or exogenous, production of ROS species typically leads to increased oxidation of NADPH, but also high flux rerouting where J ppp > J pgi [18,33]. In this physiological context, H 2 O 2 , a major source of ROS, directly interacts with and inhibits several glycolytic enzymes, notably GAPD and PFKFB3, through S-gluthationylation modifications. We therefore apply now metabolic control analysis to a context where the increase in GR-dependent oxidation of NADPH into NADP + is mediated by an increased production of H 2 O 2 , shifting the nature of parametric perturbation from k gr to k ox for which control coefficients are computed. The model incorporates now the H 2 O 2 -dependent oxidative inhibitions of GAPD and PFK1. In this scenario, control manifold equations can be derived (see Appendix C.2) to obtain a simple general expression for the NADPH control coefficient: To analyze the control properties associated with this expression, random sampling of kinetic space is performed again to plot the statistical distribution of C S,nh ox as function of key regulatory and steady-state parameters of Equation (6). The sampling procedure is supplemented by the criteria that S nh = S n to leave aside the previously characterized effect of NADPH concentration. Results show that the oxidative inhibition of glycolytic enzymes by H 2 O 2 promotes significant NADPH homeostasis, but only for negative enough level of J pgi (compare Figure 4A (6). It is to note that inhibitions of PFK1 and GAPD exhibit qualitatively the same control effect. The interplay between the feedback control of G6PD and the inhibition of lower glycolysis is depicted in Figure 4E, where the upper bound for C S,nh ox increases with both r 1 and C S, f 6p gr in independent manner. This interplay nevertheless requires a large bidirectional flux in PGI reaction as compared to the net flux (J + pgi + J − pgi J pgi ) as shown by plotting control coefficients onto the control manifold associated with Equation (6) ( Figure 4F).
To summarize, the upregulation of G6PD enzyme mediated by r 1 and the inhibition of lower glycolysis mediated by r 4/5 do not display a synergistic effect, but a strong complementary effect as these two classes of regulation promote NADPH homeostasis for, respectively, J ppp low and high compared to J pgi ( Figure 4G). Because these two classes of regulation are efficient for large J − pgi and large J + pgi , respectively, their complementary effect is empowered by high bidirectional flux (J +/− pgi > J pgi ), supporting the need for nearequilibrium PGI reaction associated with a low Gibbs free energy ∆G = J pgi ln (J + pgi /J − pgi ).  (6) and (A12) (colored arrows) recapitulating the interplay of r 1 , r 4 and r 5 and key steady-state variables on the NADPH homeostasis control coefficient.

Discussion
The present study proposes a methodological framework to investigate the mode of regulatory crosstalks in the functional response of metabolic networks. The contextdependent control of metabolic fluxes involves complex regulation patterns as exemplified in the regulation of carbon flux rerouting into oxidative branch of the PPP to meet cellular demand in NADPH [16,22]. A number of experimental and computational studies has disputed the respective or prevalent roles of diverse metabolic regulations prone to contribute to such flux rerouting [17,27,31,34,35]. The present framework aims to reconcile these studies by carefully addressing context-dependency and cooperation in the regulation pattern of NADPH homeostasis. Explicit mathematical expression of control coefficients allows us to highlight to which extent a regulatory feedback may improve, or not, a metabolic function, but more importantly how such effect depends on steady-state variables and the presence of other regulations. This provides a refined picture of complex regulation patterns that disentangles several modes of regulatory interplay ranging from synergistic effects, compensatory effects or complementary effects. Synergistic allostery usually refers to allosteric binding to separate sites of a same enzyme thereby increasing enzymatic activity in a cooperative manner [36,37]. The synergistic effect described here rather characterizes a cooperative scheme where a control coefficient associated with a given metabolic functionality is enhanced congruently by first-order effects of regulations and second-order effects where one regulation triggers concentration changes in effectors of other regulations. Such an effect falls within higher-order approaches of metabolic control analysis where, for instance, second-order control coefficients describe synergies between enzyme pairs [9,38]. The careful analysis of the cooperation between NADPH-dependent inhibition of G6PD and 6PG-dependent inhibition of PGI depicts the set of requirements for efficient synergy. For instance, increase in 6PG levels, commonly observed in response to oxidative stress [17,31], serves as a proxy to support synergistic allosteric regulation, and such an increase absolutely requires differential elasticities properties in the two main branches of the oxidative PPP controlled by G6PD and 6PGD. Specifically, 6PGD enzymes must be characterized with a lower NADPH-dependent feedback inhibition than G6PD. Although experimental measurements strongly depend on cell type and methods, comparative studies seem to indicate a larger K i value about ≈30 µM in 6PGD [39,40] compared to G6PD about ≈7 µM [25]. Rapid post-translation modifications upregulating specifically G6PD enzymes [41] could also contribute to differential upregulation of G6PD compared to 6PGD, required for 6PG increase and concomitant inhibition of PGI.
Control coefficients computed for a reference metabolic state can already identify a set of rate-limiting enzymes that are likely to be conjointly regulated supporting the requirement of distributed control [42,43]. In this scheme, the control of the regulator species itself is treated separately and can be metabolites or signaling hub proteins involved in metabolic regulation such as AMPK, AKT, or NRF2 [44][45][46]. However, an approach addressing how control coefficients depend on the reference metabolic state reveals more elaborate coordination between regulation, based on the notion of complementary regulatory efficiencies. In the case of carbon flux rerouting toward oxPPP, our results show that synergistic coregulation of PGI and G6PD is the most efficient for a metabolic network operating in a pure glycolytic mode (glyclolysis flux much larger than oxPPP flux) while the inhibition of lower glycolysis is the most efficient in the case of a PPP mode or reversed glycolytic mode. The two regulation strategies can therefore relay and compensate each other in a context-dependent manner for instance related to the severity or the source of the metabolic stress. Furthermore, such complementary range of regulation efficiency requires some particular kinetic features of the metabolic network, as PGI reactions must operate close to equilibrium. This thermodynamic requirement coincides with recent experimental results that emphasize a near-equilibrium activity of PGI enzymes in various contexts [47], including oxidative bursts [48]. This latter result motivates us to refine our framework to integrate thermodynamic constraints as already implemented in some previous studies [5,7,9].

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Kinetic Model
The metabolic network of Figure 1A describes a simplified scheme of the pentose phosphate pathway comprising n m = 9 metabolite species and n r = 12 reactions. The corresponding ordinary differential equation system ds/dt = Nv reads: The simplifying assumptions consist of pooling or neglecting some reactions as compared to a detailed kinetic model of PPP [18]: • RPI, RPE, TKT1, TKT2, and TAL are pooled to form a single reaction for the nonoxidative branch of the PPP (S7P, E4P metabolites are not included). • ALD and TPI are pooled (DHAP metabolite is not included). • GP and GRX are pooled (GSSG, GSH are not included). • Catalase reaction degrading H 2 O 2 is neglected.
We consider the conservation relation s n + s nh = 1, defining a concentration unit. Defining N + and N − as matrices collecting positive-valued stochiometric coefficients of the reactants and of the products, respectively, such that N = N − − N + , a general form of reaction rates v i reads for regulated reversible reactions: where the second term is null for irreversible reactions (i.e., K i very large) and s k is the concentration of the inhibitor associated with regulation r i . Note that v hk is a constant parameter. At steady-state concentration S, we have J i ≡ v i (S) and for reversible re- i are the first and second terms of Equation (A2). Parameter sampling is made on the parameter k i and K i randomly chosen following a uniform distribution between 10 −2 and 10 2 centered around a reference value p re f = [1, 1, 0.5, 1, 1, 2, 2, 0.5, 0.5, 1.3, 0.5, 2, 1, 2] T associated with a biologically-reasonable flux state and concentration state where:

Appendix B. Matrix Equation for Control Coefficients
For the metabolic network considered in the present study, the matrix relation C J = I + C S becomes a linear equation system where we focus on the control coefficients with respect to a given perturbation x: The perturbation corresponds either to increased NADPH consumption, (i.e., x = k gr ) or increased ROS production (i.e., x = k ox ). We useJ +/− i = J +/− i /J i in Equation (A4) only, for notation conciseness. Elasticity coefficients are given by: These elasticities associated with inhibitory interactions are negative-valued ( < 0), with minimal value of −1 for 3,4,5 . The calculations in Appendix C may use the asymptotic property that lim S nh →0 1,2 = S nh (1 + r 1,2 ).
For a given set of parameters p * , the steady state equation N v(S(p * ), p * ) = 0 has a unique solution {S * , J * } = f (p * ). Associated to a particular steady-state, the above equation thus gives a unique solution for control coefficient C J,S = f (p * ) which are computed during the parameter sampling procedure.

Appendix C. Computation of Control Manifolds
Appendix C.1. Control Analysis of Regulatory Crosstalk r 1,2,3 A first setting of the problem focuses on the interplay between the regulation of oxPPP enzymes and PGI (r 1,2,3 in Figure 3). We consider a perturbation of the oxidation rate of NADPH where x = gr, π gr = [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0] T and a gr = 0 in Equation (A4). We define J ppp as the steady-state flux through the oxidative branch of the PPP, which relates to some other steady-state fluxes as: Combining Equation (A8), a more general control expression is obtained: Assuming a pure glycolytic mode J pgi = J + pgi = J hk , we obtain a set of asymptotic relations: max C J,ppp gr • All those asymptotic relations are proportional to S nh , justifying the use of the normalized control coefficient C J,ppp gr /S nh (Figure 3). • Equation (A11a) shows that a PPP flux control driven by N ADP + cofactor binding to G6PD (no regulation) has an upper bound of S nh . • Equation (A11b) expresses that r 1 promotes PPP flux control (i) independently on r 2 , (ii) especially for small S nh . • Equation (A11c) defines the complex nonlinear interplay between contributions of r 1 , r 2 , and r 3 . • Equation (A11d) shows indeed that r 3 alone, even very large, cannot increase PPP flux control. • Equations (A11e) and (A11f) indicate a synergistic cooperation between r 1 and r 3 (product term) for both small or large values of r 3 .
Beyond the assumption of a pure glycolytic mode, Equation (A10) contains the terms J ppp in the denominator and J pgi in factor of 3 , which means that flux control decreases with increasing value of J ppp and that the effect of 2 cancels for J pgi = 0 (i.e., pure oxidative PPP mode).
Appendix C.2. Control Analysis of Regulatory Crosstalk r 1,4,5 An alternative setting of the problem focuses on the interplay between the regulations of oxPPP enzymes and PFK1 or GAPD (r 1,4,5 in Figure 4). The perturbation of NADPH homeostasis is now due to an increased production rate of H 2 O 2 or ROS (i.e., k ox ) where π ox = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1] T and a ox = 1 in Equation (A4). For such perturbation scheme, the relevant set of control equations becomes: where one obtains after combining Equations (A12a)-(A12d): which finally simplifies for 3 = 0 to: (A14) NADPH homeostasis is related to smallest |C S,nh ox | values as possible and can be promoted by several mechanisms: • Positive values of C S, f 6p ox which is increased by regulation r 4 following Equation (A12f), but also regulation r 5 . • The effect of those regulations is enhanced by large directional PGI flux J − pgi from F6P to G6P. • The effect of r 1 is amplified by large directional PGI flux J + pgi from G6P to F6P. • The two items above indicate that efficient regulation of NADPH homeostasis by r 1,4,5 requires high values of both J +,− pgi with an upper bound given by :