Differential suitability of reactive oxygen species and the role of glutathione in regulating paradoxical behavior in gliomas: A mathematical perspective

Manipulative strategies of ROS in cancer are often exhibited as changes in the redox and thiol ratio of the cells. Cellular responses to oxidative insults are generated in response to these changes which are triggered due to the rerouting of the metabolic framework to maintain survival under stress. However, mechanisms of these metabolic re-routing are not clearly understood and remained debatable. In the present work, we have designed a context-based dynamic metabolic model to establish that the coordinated functioning of glutathione peroxidase (GTHP), glutathione oxidoreductase (GTHO) and NADPH oxidase (NOX) is crucial in determining cancerous transformation, specifically in gliomas. Further, we propose that the puzzling duality of ROS (represented by changes in h2o2 in the present model) in exhibiting varying cellular fates can be determined by considering simultaneous changes in nadph/nadp+ and gsh/gssg that occur during the reprogramming of metabolic reactions. This will be helpful in determining the pro-apoptotic or anti-apoptotic fate of gliomas and can be useful in designing effective pro-oxidant and/or anti-oxidant therapeutic approaches against gliomas.


Rate Equations in Cleland Nomenclature
Kinetic equations for all the metabolic reactions were written following the Cleland nomenclature for mono-substrate, bi-substrate, and ter-substrate reactions. Kinetic rate equations for the enzymes have been previously derived mechanistically and reported in various literature. Initial kinetics with first-order reaction rates of the enzymatic reactions have been considered for writing the equations. Since initial rates are usually measured in the presence of substrates only, there are no product terms present in the equations. Values for all the parameters have been listed in Table S1 and initial values of all the variables have been provided in section 4 of the supporting information.

Central Carbon Metabolism
Glycolysis: Kinetic rate equations for the enzymes belonging to the first half of glycolysis, from where it branches to the glycine-serine pathway have been considered in the model. The kinetics is mostly represented by uni-uni and bi-bi mechanisms. The equation forms for most of the reactions are similar except for GAPDH and PGK. Competitive inhibition of GAPDH by h 2 o 2 has been considered. The rate kinetics of PGK follows a partial rapid equilibrium random bi-bi steady-state kinetics [1] where the binding of 1,3bpg with the enzyme is transient and has not been considered in the denominator. The assumption here is that the final products of the pathway, i.e., fructose and glyceraldehyde 3-phosphate re-enter the glycolytic pathway. The whole pathway has been reduced to a single equation and the parameter values for the equation have been determined by parameter estimation technique to represent the real biological situation. Competitive inhibition of the pathway by atp has been considered. Eq.9:

Differential Equations:
All the differential equations with metabolites as the variables have been formulated by considering the involvement of the metabolites in different enzymatic reactions. A few natural productions and decay rates have been considered to take into account the formation of metabolites from other sources and their utilization in other cellular processes respectively. Values to all the parameters considered in the differential equations have been tabulated in Table  S1.
gly e d gly l v d gly dt

Positivity and Boundedness
Positivity:   R . In this region, the usual existence, uniqueness and continuation results hold for the system. From our numerical simulations, we have observed the existence of positive solutions for the parameter values presented in Table S1. The steady-state attained by all the variables of the model is a part of this positive solution space.

Parameter values and reported range of initial values of variables:
All the parameters have been defined and the values have been provided along with the reference to the literature from which they have been obtained. A set of parameters has been estimated using the Delayed Rejection Adaptive Metropolis (DRAM) algorithm of Markov Chain Monte Carlo (MCMC) Toolbox in MATLAB 2017a. These parameters have been referenced as "Estimated". Another set of parameters has been referenced as "Expected". These parameters have been estimated by varying the parameters within the biologically feasible ranges as reported in various literature to determine their expected value to calibrate the model with experimental observations.  [11] 60.
GCL glut k Rate constant for association of Glutamate with Glutamyl-cysteine Ligase 0.46 mM [11] 61.
GS gly k Rate constant for association of Glycine with Glutathione synthase 1.37 mM [13] 68.
GTHO nadph k Rate constant for association of NADPH with Glutathione Oxidoreductase 0.063 mM [21] 87.

Parameter estimation:
Parameter estimation of 23 unknown parameters was done using the Delayed Rejection Adaptive Metropolis (DRAM) algorithm of Markov Chain Monte Carlo (MCMC) Toolbox in MATLAB 2017a. The parameter distribution plot of the 23 parameters has been shown in Figure S1A and the trace plot till 5 lakh chains has been shown in Figure S1B. The model was re-evaluated for available data for glutamate exchange in astrocytes [45] to estimate the value of parameters GLUTEX m V and GLUTEX m k . The predictive plot for the same has been provided in Figure S2A and the distribution plot and trace plot of estimated parameters have been shown in Figure S2B and S2C respectively.

Sensitivity Analysis:
Extended Fourier Amplified Sensitivity Test (eFAST) algorithm was used for identifying sensitive parameters to the system. The sensitivity analysis was carried out using the whole set of parameters under normal glial, hypoxic and glioma scenarios respectively. The sensitivity indices of the parameters were captured for multiple time points (1hr, 2hrs, 4hrs, 8hrs, 16hrs, 20hrs, 24hrs, 32hrs, 50hrs, 100hrs, 150hrs, and 200hrs) and the difference in Si (if any) was taken into account.   Combinatorial variations of sensitive parameters reported in the glioma scenario G11 specified in Table 1 of the main text were made. The values of the parameters were varied between a wide range (0.0001 to 100 units) and changes were observed. Here we have reported the temporal variation in the h 2 o 2 , gsh/gssg and nadph/nadp + profiles for a particular value of the parameters. The values of the parameters have been specified in Table S2 and the temporal plots have been shown in Figure S6.