A Dynamic Metabolic Flux Analysis of Myeloid-Derived Suppressor Cells Confirms Immunosuppression-Related Metabolic Plasticity

Recent years have witnessed an increasing interest at understanding the role of myeloid-derived suppressor cells (MDSCs) in cancer-induced immunosuppression, with efforts to inhibit their maturation and/or their activity. We have thus modelled MDSCs central carbon metabolism and bioenergetics dynamic, calibrating the model using experimental data on in vitro matured mice bone marrow cells into MDSCs. The model was then used to probe the cells metabolic state and dynamics, performing a dynamic metabolic flux analysis (dMFA) study. Indeed, MDSCs maturation correlates with a high glycolytic flux contributing to up to 95% of the global ATP turnover rate, while most of the glucose-derived carbon enters the TCA cycle. Model simulations also reveal that pentose phosphate pathway and oxidative phosphorylation activities were kept at minimal levels to ensure NADPH production and anabolic precursors synthesis. Surprisingly, MDSCs immunosuppressive activity, i.e. L-arginine uptake, metabolism and endogenous synthesis, only consumes sparse quantities of energy-rich nucleotides (ATP and NADPH). Therefore, model simulations suggest that MDSCs exhibit a heterogeous metabolic profile similar to tumour cells. This behavior is probably an indirect immunosuppressive mechanism where MDSCs reduce the availability of carbon sources in the tumour periphery microenvironment, which could explain the dysfuntion and death of immune effector cells.

Dynamic metabolic flux analysis (dMFA) have shown being useful, as a quite unique tool, to describe extracellular and intracellular nutrients and metabolites concentration as well as fluxes transient behavior in mammalian cells [9][10][11] . However, such approach, as well as other metabolic modelling approaches, have not been applied yet to study the immunosuppression phenomenon, to the best of our knowledge. We thus propose here to apply a modelling approach [9][10][11] previously validated on Chinese Hamster Ovary (CHO) cells [9][10][11] , on the maturation process of bone marrow cells into MDSCs, and thus perform dMFA to unreaveal metabolic traits of MDSCs' maturation process.
The dynamic metabolic model combines the central carbon metabolism (CCM) and cell energetics, redox states and metabolic regulation mechanisms, and includes MDSCs immunosuppressive machinery (ARG1 and iNOS enzymes) and L-arginine endogenous synthesis reactions. The model was calibrated on experimental data of extra-and intracellular metabolites generated in a previous work 8 . We are thus first presenting a descriptive model as well as evaluating its simulation capacity. Sensitive kinetic parameters and their confidence intervals were determined to identify the key kinetic parameters and biochemical reactions to be considered as potential targets for modulating MDSCs maturation and immunosuppressive activity. Then, a dMFA study was performed through model simulations for evaluating MDSCs metabolic performance.
Model simulations confirmed the hypothesis that MDSCs maturation correlates with a high glycolytic flux sustaining ATP turnover and also supporting a high TCA cycle activity, with low but not suppressed pentose phosphate pathway and oxidative phosphorylation activities generating NADPH and anabolic precursors. Interestingly, MDSCs immunosuppressive machinery reveals not being energetically costly, while the overall nutritional profile and metabolic flux studies showed that MDSCs mimic tumour cells' metabolic plasticity 12,13 .

Materials and Methods
Cell Culture. Experimental data used herein to verify the accuracy of the metabolic model were taken from our previous work 8 where experimental settings are well detailed.
Briefly, BM cells were extracted from 6-to 8-weeks old C57BL/6 mice (Charles River, Quebec, Canada). Animal experimentations were performed in accordance with the Canadian Council on Animal Care guidelines and the protocol was approved by the Université de Montréal's ethical committee. Mice were kept under specific pathogen-free conditions prior to euthanasia by CO 2 .
It should be noted that the term "BM-derived MDSCs" is used when referring to experimental and simulation results on BM cells derived into MDSCs. When referring to literature, the more general term MDSCs is used because it is the terminology normally accepted for these cells. Metabolites Measurements. Every 24 hours, extracellular nutrients and metabolites in medium were measured using a dual-channel immobilized oxidase enzyme biochemistry analyzer (2700 SELECT, YSI Life Sciences, USA) and metabolites were extracted using cold methanol and sonication on ice method modified from Kimball et al. 8,15 .
Nucleotides concentrations were measured using a 1290 UPLC system coupled to a 6460 triple quadruple mass spectrometer. A Symmetry C18 column (Waters, Canada) was employed for the separation. Mobile phases consisted of 10 mM ammonium acetate with 15 mM DMHA (N,N-dimethylhexylanine) at pH 7.0 and 50%/50% (v/v) acetonitrile, 20 mM NH 4 OAc at pH 7.0. Organic acids concentrations were assessed using the       Model Development and Calibration. The construction steps while building the structure of the metabolic model presented in this study is based on previous models describing the central carbon metabolism of plant 16,17 and CHO cells 10,11,18 where it showed satisfactory results simulating experimental data. The global procedure for model development is described as Supplementary material. Briefly, we have here developed a model integrating metabolic pathways known to be expressed in MDSCs, by applying a modelling approach thus previously validated with eukaryotic cells. Basically, the network of biochemical reactions was first built including central carbon metabolism as well as biochemical reactions known to be of major importance in MDSCs. Metabolic flux kinetics and regulation mechanisms were then described, based on literature on MDSCs when available, or on previous works on animal cells. The model was then describing known metabolic network in MDSCs as well as it was anchored specifically onto MDSCs by performing model structure calibration and parameters value identification thanks to experimental in vitro metabolomic data with these cells.
The model metabolic network (Fig. 1) (see nomenclature in Table S1) represents the main biochemical pathways in mammalian cells including glycolysis, glutaminolysis, pentose phosphate pathway, TCA cycle, cell energetics and oxidative phosphorylation (i.e. respiration) (Fig. 1A). In addition, the urea cycle (Fig. 1B) was added in order to model the reactions involving the L-arginine metabolism, which is specifically active in the MDSCs immunosuppressive phenomenon 3 . Furthermore, the catabolic pathways related to the amino acids metabolism (notably alanine, asparagine and aspartate) were also considered since these are known to contribute as carbon source. Therefore, the model describes major inputs, such as major nutrients, and outputs, such as cell growth and the management of cell metabolites from the immunosuppressive phenomenon. This biosystem network is thus highly simplified but it considers the major set of interlinked biochemical reactions allowing to describe MDSCs behaviour.
Moreover, the modelling strategy is based on the following considerations/assumptions: • Flux kinetics consider the effect of substrates as well as co-factors (i.e. energetic and redox nucleotides).
• Except for the amino acids alanine, asparagine, aspartate, glutamine and glutamate, all other extracellular metabolite concentrations are considered higher that the affinity constants of associated cell membrane transporters 19 . Therefore, intra-and extracellular metabolites that are consumed or secreted were taken as being at the same concentration across the cell membrane; a reliable assumption we have demonstrated in past works dealing with similar models 10,11,[16][17][18] . • The reactions of the adenylate kinase and ATPase enzymes are used to balance the time evolution of AMP, ADP and ATP energetic shuttles, with the synthesis of AMP being described. • The model assumes non-constant energetic metabolites concentration: NAD + , NADH and NADP + , NADPH.
NAD + synthesis is thus described, from which NADP is synthesized.  Table 3. Metabolites initial concentration.
• Enzymatic activity of the NADPH oxidase and the mitochondrial proton leak phenomenon are employed to balance the time evolution of NADP + /NADPH and NAD + /NADH, respectively. • The stoichiometry of the reaction describing cell division and growth was based on average cell composition for CHO cells, since this information for MDSCs was not available 9 .
The final (i.e. after performing sensitivity analysis -model parameter values optimization cycle) dynamic model includes 44 reactions (Table 1) and mass balances (Table 2), with a total of 136 kinetic parameters ( Table 2). The stoichiometric coefficients of the respective biosynthetic equations listed in Table 1 were taken from literature 20 . The fluxes are modeled by multiplicative Michaelis-Menten type equations and a factor of non-competitive inhibition is added when metabolites show inhibitory effects. Nutrients and metabolites as well as cell initial concentration (i.e. at t = 0) are presented in Table 3.
Using the parameter values determined in a previous work on CHO cells 11 (and from references therein), a first sensitivity analysis was performed on the resulting model, aiming to identify the most critical (i.e. sensitive) parameters. In brief (see the Supplementary material for details), because of the large number of parameter values to be determined, an iterative procedure including manual trials on sub-groups of parameters and a model simulation global error minimization step, allowed identifying optimal model parameter values. Sensitivity analysis of model parameters on the model acuity to simulate experimental reality was performed as follow. Values of model parameters were changed by ±5 or 10%, once at a time, from their optimal value, and the normalized sum-squared differences were calculated as described below. A model sensitivity analysis on initial conditions has also been performed to evaluate the effect of measurement errors. X mea and X sim are respectively the experimental data and simulated values for each state variable m and time k, and the weight is the inverse of the variance of the experimental data for each state variable, − var m 1 .  Finally, the 95% confidence intervals on the estimated model parameters were computed, using the Matlab functions "nlinfit.m", "nlparci.m" and "nlpredci.m". Following the parameter estimation step ("nlinfit.m"), the confidence intervals were computed for the estimated model parameters ("nlparci.m"). The 30 most sensitive parameters (Table 4) have been used in the computation of the confidence intervals.

Results and Discussion
The Kinetic-Metabolic Model Simulates MDSCs Metabolic-Time Profile. The resulting metabolic model has then been used to simulate metabolites concentration with time, supported (Fig. 2) or not (Fig. S2) with experimental data. The time-evolution of metabolites concentration were thus either compared to experimental data (when available) and to range of values found in literature, and thus confirmed the biologic relevance of the set of parameters as well as of the metabolic model structure.
The sensitivity analysis (Fig. 3A) procedure allowed to rank the parameters from their decreasing influence, and to remove parameters that were not contributing to model sensitivity from further optimization cycles, keeping them at their optimal values. Only the sensitivity analysis results for the final model (i.e. no more improvement of the simulation error) are shown here. Interestingly, most of the parameters (102 on 131) reveal to be non-sensitive with a WSSRES value lower than 0.1. This lack of sensitivity may partially come from the large number of parameters to be optimized, the low number of experimental data, and the important error bars. These non-sensitive parameters are within biologically relevant values and describe existing active pathways and enzymatic reactions, but they may require an expanded experimental space out of the current one to be solicited. Moreover and of interest, we have previously showed that the sensitive parameters are the ones needing to be changed while describing a modified genetic expression 10,11 . Taking all of the above, the model structure was not reduced. Specifically, the model reveals to be primarily sensitive to parameters of glycolysis and energetic reactions, partially to TCA cycle and pentose phosphate pathway, and to a lesser extent to urea cycle and amino acids catabolism. The maximum specific glucose uptake rate (ν maxHK ) and the maximal ATPase (ν maxATPase ) reaction rate as ) showed to strongly affect simulation error. This result is particularly interesting since cell energetics plays a major role on the regulation of central carbon metabolism fluxes 18 . In complement, we also observed that the model is less sensitive to errors on initial conditions (Table 3) (Fig. 3B), only showing here the ones with the highest effect on the global error. This result is of interest since this type of dynamic model only requires initial conditions, together with kinetic parameters, to simulate a cell population behavior with time. Indeed, it is interesting to note that the high simplification level of the proposed model is not resulting in high level of simultation errors when challenged from either parameters value or initial concentration conditions. Expanding the amount of biochemical reactions would require, in a further work, expanding the amount of quantified metabolites to support the according increase of model parameters. The model has thus been used to perform a dynamic metabolic flux analysis (dMFA), analysing the results as new data per se as well as to further validate model simulations validity.

MDSCs Exhibit Warburg Effect During their Maturation Process. Dynamic metabolic flux analysis
confirms MDSCs mimick tumour cells metabolic profile. Indeed, model simulation of experimental data confirm that BM-derived MDSCs exhibit a high glycolytic activity level as revealed by the continuous increase of the flux through pyruvate kinase (Fig. 4A). Moreover, up to 95% of the ATP generated in MDSCs are glycolysis-dependent (Fig. 4B), with thus an expected reduced oxygen consumption rate. This result correlates with previous findings where BM-derived MDSCs were shown to decrease their oxygen consumption, and consequently the oxidative phosphorylation -dependent ATP production, by 60% during their maturation process 8 . This behavior is typical to tumor cells. In fact, tumour cells are known to be metabolically heterogeneous 12, 13 with a high glycolytic metabolism and a reduced, but not completely suppressed, oxidative phosphorylation 21 . A high aerobic metabolism was also shown to support highly proliferative tumors cells whereas a less efficient electron transfer chain in mitochondria is associated to ROS production 22 . This result is of utmost interest since despite the abundance of oxygen in in vitro culture condition, a glycolytic metabolism is favored in BM-derived MDSCs, a behaviour reported in some tumour cells. This metabolic profile may represent another mechanism of immunosuppression. In fact, in the tumour microenvironment, anti-tumour immune effector cells will compete with tumour cells and MDSCs that both exhibit high glucose uptake rates. Whereas, immune cells do not have any metabolic elasticity to acclimatize to low oxygen tension and limited glucose availability, and these conditions may trigger immune cell dysfunction and death, which can indirectly lead to tumour escape and progression. In a previous work on BM-derived MDSCs and also on a MSC-1 immortalized cell line 8 , we observed several similarities  between MDSCs and tumor cells metabolic behavior, such as a glycolytic metabolism, high glucose and glutamine uptake rates, reduced oxygen consumption rate but a high TCA cycle activity. It is known that MDSCs are getting fully activated and matured in the tumor microenvironment by means of tumor-derived factors (TDFs such as GM-CSF, TGF, IL-6, etc). These TDFs indirectly induce a similar metabolic behavior to tumor cells and this is further confirmed in this work from the dynamic metabolic flux analysis (dMFA) performed using our model. ). (C) NADPH-to-NADP ratio. (D) Contribution of L-glutamine to TCA cycle intermediates replenishment (

. (B) Glycolysis contribution to ATP production
). All results obtained from the same model simulation than that of Fig. 2. Most of the Glucose-Derived Carbon Enters the TCA Cycle. A dMFA study at the pyruvate branch point revealed that pyruvate is mainly produced through the pyruvate kinase enzyme, with a flux of ~700-fold than that through the malic enzyme. This result suggests that pyruvate is principally derived from glucose than recirculated from the TCA cycle (Fig. 4C). Likewise, the comparison of the fluxes allowing the entry of glucose-derived carbon into the TCA cycle showed that the flux through pyruvate dehydrogenase was higher (up to 50-fold), and kept increasing with time, than the flux through pyruvate carboxylase (Fig. 4C). Moreover, the percentage of glucose-derived carbon entering the TCA cycle was increasing (from 50 to 80%) during the BM-derived MDSC maturation process (Fig. 4D). These anaplerotic/cataplerotic reactions are known to be important to the replenishment of TCA cycle intermediates 23 . The high glycolysis-TCA cycle activity is probably associated to the acquisition of the immunosuppressive phenotype, and the maintenance of inflammatory environment at the tumour edge, by supporting enzyme and surface markers expression, cytokine/chemokine synthesis, etc. In addition, some of the TCA cycle intermediates have roles other than being a carbon source. For instance, the accumulation of fumarate and succinate stabilizes hypoxia induction factor 1α 22, 24 that was shown to be expressed in MDSCs to adapt to low oxygen tension in tumour 25 . Since BM-derived MDSCs exhibit a low growth rate during the 96 hours of the study, a longer culture period or the analysis of flux distribution in in vivo generated MDSCs is required to confirm if this behavior, i.e. high TCA cycle activity, is associated to proliferative purposes as previously reported in cancer cells 22 . Moreover, the analysis of the G6P branch point shows that only 20% of the glucose-derived carbon enters the pentose phosphate pathway, since the ratio of the flux through G6PDH to the one through HK was kept constant at 20% (Fig. 5A). A total of 30% of the carbon that enters the pentose phosphate pathway is recirculated to glycolysis (Fig. 5B). This metabolic behavior occurs when both NADPH and ATP are needed, but only sparsely ribose-5-phosphate, which is coherent with the low growth rate observed 8 . Glycolytic carbons are then shunted into the oxidative phase of the pentose phosphate pathway, and consequently in the non-oxidative phase leading to re-entering glycolysis 26 . This was confirmed by the NADPH-to-NADP ratio (Fig. 5C) which increased (up to 9) rapidly for the last 30 hours, where BM-derived MDSCs were fully mature and active. Furthermore, as we discussed above, glycolysis was mainly responsible for ATP production, which is thus confirming then that MDSCs have the ability to modulate their central carbon metabolism to ensure the proper bioenergetics state required for their maturation and activity.

MDSCs Immunosuppression Machinery is not Energetically Costly.
Contrarily to glucose-derived carbon, which is mainly entering into the TCA cycle to sustain its high activity, L-glutamine (GLN in the model) only contributes sparsely to the replenishment of TCA cycle intermediates (Fig. 5D). However, the fate of GLN-derived carbon in MDSCs is ambiguous and requires further investigation using 13 C-labelled GLN, but our results clearly show it is likely used in the endogenous synthesis of L-arginine to support cells immunosuppressive potential 27 . As BM-derived MDSCs exhibit a high GLN consumption rate, immune effector cells will also compete with MDSCs in addition to tumour cells for GLN, a phenomenon that is crucial for their proliferation.
Interestingly, the dMFA study showed that the L-arginine metabolism (including its uptake, its metabolism by the inducible Nitric Oxide Synthase and Arginase 1 and its endogenous synthesis) do not use considerable quantities of ATP and NADPH. In fact, the percentage of NADPH used by the L-arginine metabolism is decreasing from 5 to 1% (Fig. 6A). This decrease is associated to a higher flux through NADPH oxidase during BM-derived MDSCs maturation process. However, NADPH oxidase co-produces NADP + and superoxide (O 2 − ), the latter being a precursor for the generation of immunosuppressive species ROS and RNOS 28 . A similar decreasing trend from 6 to 1% was simulated for the percentage of ATP consumed by the L-arginine metabolism (Fig. 6B). These findings support the hypothesis that L-arginine metabolism is not energetically costly but further investigation is still required to verify the energetic needs of complementary immunosuppressive mechanisms that are not modelled in this work, such as tryptophan metabolism, expression of messenger molecules, etc. Interestingly, there are only 30 sensitive parameters (Table 4) on the 136 parameters of the model, which cover all metabolic sub-networks, and especially including those that define specific functions of MDSCs.

Conclusion
This work on model simulation of MDSCs maturation, used a wide set of experimental data for extra-and intracellular metabolites concentration to develop a descriptive and predictive dynamic metabolic model. The model then enabled performing a dMFA study, which results allowed highlighting the carbon distribution dynamics during MDSCs maturation and also in fully mature and active cells. Mainly, MDSCs revealed to be dependent on the glycolysis and the glucose-derived carbon to ensure high central carbon metabolism activity and sustained production of ATP. Pentose phosphate pathway and oxidative phosphorylation showed being active at the minimal level allowing to replenish the NADPH pool and anabolic precursors. This metabolic behaviour was not justified by specific requirements of MDSCs immunosuppression machinery neither at the carbon intermediate or energy level. However, this metabolic profile indirectly favors tumour invasiveness and shall be considered as a target for anti-cancer immunotherapy. In this work, the interpretation of carbon flux distribution simulated by the model leads to the identification of the specific metabolic signatures of mature MDSCs that are directly associated to the immunosuppression phenotype. This in silico platform thus allows prospecting MDSCs dynamic behaviour under stimulation or inhibition conditions, as well as it represents a complementary tool to existing tools for studying the immunosuppression phenomenon 29 . Cell metabolomic studies provide new and complementary insights enabling to better understand the immunosuppression phenomenon, and such dynamic model may be useful to the identification of novel biomarkers as well as to define therapeutic strategies to modulate this phenomenon. The validation of the predictive capacity of the model will, however, be possible with the accumulation of large experimental data sets, as well as further