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.

Introduction

The invasiveness of conventional anti-cancer therapies appeals to the establishment of novel therapeutic strategies. Specifically, the immunotherapy approach is gaining in interest as it showed promising results in the treatment of cancer of various types[1]. One of the mechanisms of immunotherapy specifically targets the inhibition of myeloid-derived suppressor cells (MDSCs) maturation and activity. Briefly, MDSCs are a heterogenous cell population of macrophages, dendritic cells, granulocytes, to name only the most abundant cell types[2]. These cells metabolize the semi-essential amino acid L-arginine through the enzymatic activities of arginase 1 and inducible nitric oxide synthase (ARG1 and iNOS, respectively)[3]. Both the sparse availability of L-arginine and the accumulation of nitric oxide derivatives (e.g. reactive nitrogen oxide species, RNOS) trigger the suppression of specific anti-tumour immune response[4]. Despite the considerable progress achieved in the comprehension of MDSCs biology that has resulted in multiple therapeutic approaches proposed to overcome MDSCs-mediated tumor escape, the metabolic events supporting MDSCs immunosuppressive potential are still ambiguous and poorly investigated.

Previous work on the activation of dendritic cells with lipopolysaccharides (LPS, known to activate iNOS enzyme) showed that cells undergo a PI3K/Akt-dependent metabolic transition from oxidative phosphorylation to aerobic glycolysis[5]. Likewise, LPS-activated M1 macrophages display enhanced glycolytic metabolism and reduced mitochondrial activity[6]. Therefore, we concentrated our efforts to study the central carbon metabolism and bioenergetics of MDSCs, investigating their specific requirements to acquire an immunosuppressive phenotype[7, 8]. The analysis of intracelullar metabolites concentration suggested that MDSCs exhibit a glycolytic metabolism and a high TCA cycle activity providing the immunosuppressive mechanisms with energy-rich nucleotides and carbon intermediates.

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-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-11] previously validated on Chinese Hamster Ovary (CHO) cells[9-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 CO2.

BM cells were cultured in RPMI1640 medium (Sigma, Ontario, Canada) supplemented with 10% (v/v) irradiated fetal bovine serum (Cedarlane, Onatrio, Canada), 1 mM Sodium Pyruvate (Sigma), 50 µM β-Mercaptoethanol (Sigma), 100 U/mL Penicillin, 150 U/mL Streptomycin (Cedarlane) in a 5% CO2 and 37 °C incubator. The maturation of BM-derived MDSCs was performed by culturing BM cells for 96 hours in the supplemented medium mentioned above in the presence of 40 ng/mL of GM-CSF (granulocyte macrophage - colony stimulating factor) and 40 ng/mL IL-6 (Interleukin-6) as shown by Marigo et al.[14].

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 NH4OAc at pH 7.0. Organic acids concentrations were assessed using the above-mentioned UPLC-MS/MS system with a Hypercarb column (Thermo Fisher, Ontario, Canada). Mobiles phases were 20 mM amoonium acetate at pH 7.5 and 10% (v/v) methanol in water.

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.

Figure 1

Description of the metabolic network of the model. (A) Central carbon metabolism and bioenergetics. (B) Urea cycle and amino acid catabolism. Brown: extracellular, black: intracellular. Enzymes are indicated in all reactions. All reactions are explicitely described in Table 1, including the stoichiometry.

Table 1

Reactions of the metabolic network.

No.	Enzyme	Flux	Reaction
1	Hexokinase	VHK	EGLC + ATP → G6P + ADP
2	Phosphoglucose isomerase	VPGI	G6P ↔ F6P
3	Glucose-6-phosphate dehydrogenase	VG6PDH	G6P + 2 NADP+  → R5P + CO2 + 2 NADPH
4	Ribulose-5-phosphate epimerase	VEP	R5P →  × 5 P
5	Transketolase	VTK	R5P + 2 × 5 P → 2 F6P + GAP
6	Phosphofructokinase	VPFK	F6P + ATP → 2 GAP + ADP
7	Phosphoglycerate kinase	VPGK	GAP + ADP + NAD+  + Pi → PEP + ATP + NADH
8	Pyruvate kinase	VPK	PEP + ADP → PYR + ATP
9	Lactate dehydrogenase	VLDH	PYR + NADH ↔ ELAC + NAD+
10	Pyruvate dehydrogenase	VPDH	PYR + CoA + NAD+  → ACCoA + CO2 + NADH
11	Citrate synthase	VCS	ACCoA + OAA → CIT + CoA
12	Aconitase/isocitrate dehydrogenase	VCITS	CIT + NAD+  → AKG + CO2 + NADH
13	Alpha ketoglutarate dehydrogenase	VAKGDH	AKG + CoA + NAD+  → SCoA + CO2 + NADH
14	Succinyl coenzyme A synthetase	VSCOAS	SCoA + ADP + Pi → SUC + CoA + ATP
15	Succinate dehydrogenase	VSDH	SUC + 2/3 NAD+  → FUM + 2/3 NADH
16	Fumarase	VFUM	FUM → MAL
17	Malate dehydrogenase	VMLD	MAL + NAD+  → OAA + NADH
18	Malic enzyme	VME	MAL + NADP+  → PYR + CO2 + NADPH
19	Pyruvate carboxylase	VPC	PYR + CO2 + ATP → OAA + ADP + Pi
20	Alanine aminotranferase	VAlaTA	PYR + GLU → AKG + ALA
21	Glutamine transport	VGlnT	EGLN + ATP ↔ GLN + ADP
22	Glutamate transport	VGluT	GLU + ADP + Pi → EGLU + ATP
23	Glutamine synthetase	VGLNase	GLN ↔  → GLU + ENH4
24	Glutamate dehydrogenase	VGLDH	GLU + NAD+  → AKG + NADH + ENH4
25	Aspartate aminotransferase	VASTA	GLU + OAA → AKG + ASP
26	Aspartate transporter	VAspT	EASP + ATP ↔ ASP + ADP
27	Asparagine transporter	VAsnT	EASN + ATP ↔ ASN + ADP
28	Asparagine synthetase	VAsnS	ASP + ATP → ASN + AMP + PPi
29	Asparaginase	VAsnase	ASN → ASP + ENH4
30	Arginine transporter	VArgT	EARG + ATP → ARG + ADP
31	Argininosuccinate lyase	VASL	AS → ARG + FUM
32	Arginase 1	VArg1	ARG + H2O → ORN + Urea
33	Ornithine transcarbamylase	VOTC	ORN + CP → CTR + Pi
34	Nitric oxide synthase	ViNOS	ARG + 3/2 NADP+  → CTR + 3/2 NADPH
35	Argininosuccinate synthase	VASS	CTR + ASP + ATP → AS + AMP + PPi
36	Carbamoyl phosphate synthetase	VCPS	CO2 + NH4 + 2 ATP → CP + 2 ADP + Pi
37	Nucleotide synthesis	VPrAMP	R5P → AMP
38	Creatine kinase	VCK	PCr + ADP ↔ Cr + ATP
39	NADPH oxidase	VNADPHox	NADPH + 2 O2 → NADP+  + 2 O2 − + H+
40	Proton leak	Vleak	O2 + 2 NADH → 2NAD+  + 2 H2O
41	Respiration	VResp	O2 + 2 NADH + 4 ADP + 4 Pi → 4 ATP + 2 NAD+  + 2 H2O
42	ATPase	VATPase	ATP → ADP + Pi
43	Adenylate kinase	VAK	ATP + AMP → 2 ADP
44	NAD synthesis	VNAT	R5P +  × 5 P + EGLN + 2 ATP ↔ GLU + NAD + AMP
45	NADP synthesis	VNHG	NAD + ATP ↔ NADP + ADP
46	Cell growth rate	Vgrowth	α G6P + β R5P + γ EGln + δ ATP → X + δ ADP α = 3.18 × 10−7, β = 1.50 × 10−8, γ = 5.03 × 10−9, δ = 3.15 × 10−8

N.B. H+, O2, CO2, Pi, PPi and H2O were not considered in mass balances.

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–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.

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.

Table 2

Equations of kinetic fluxes and parameter values.

Number	Equation	Parameter values
1	V_HK=V_max,HKEGLCEGLC+K_1,EGLCATPATP+K_1,ATP	Vmax,HK = 5.08 × 10−4 K1,ATP = 1.08 × 10−5 K1,EGLC = 10−12
2	V_PGI=Vmax,PGIfG6PG6P+K_2,G6P−Vmax,PGIrF6PF6P+K_2,F6P	Vf max,PGI = 3.53 × 10−4 Vr max,PGI = 9.47 × 10−7 K2,F6P = 1.93 × 10−5 K2,G6P = 1.04 × 10−6
3	V_G6PDH=V_max,G6PDHG6PG6P+K_3,G6PNADPNADP+K_3,NADPATPATP+K_3,ATP	Vmax,G6PDH = 7.62 × 10−5 K3,ATP = 1.03 × 10−12 K3,G6P = 1.32 × 10−6 K3,NADP = 10−12
4	V_EP=V_max,EPR5PR5P+K_4,R5P	Vmax,EP = 6.74 × 10−5 K4,R5P = 10−6
5	V_TK=V_max,TKX5PX5P+K_5,X5PR5PR5P+K_5,R5P	Vmax,TK=3.37 × 10−5 K5,R5P = 10−6 K5,X5P = 10−6
6	V_PFK=V_max,PFKF6PF6P+K_6,F6PATPATP+K_6,ATP	Vmax,PFK = 3.42 × 10−4 K6,ATP = 1.03 × 10−12 K6,F6P = 2.83 × 10−7
7	V_PGK=V_max,PGKGAPGAP+K_7,GAPNADNAD+K_7,NADADPADP+K_7,ADP	Vmax,PGK = 7.35 × 10−3 K7,ADP = 4.75 × 10−10 K7,GAP = 10−6 K7,NAD = 10−8
8	V_PK=V_max,PKPEPPEP+K_8,PEPADPADP+K_8,ADP	Vmax,PK = 9.79 × 10−4 K8,ADP = 4.75 × 10−10 K8,PEP = 4.20 × 10−7
9	V_LDH=V_max,LDHPYRPYR+K_9,PYRNADHNADH+K_9,NADH	Vmax,LDH = 9.99 × 10−4 K9,NADH = 9.86 × 10−6 K9,PYR = 1.65 × 10−12
10	V_PDH=V_max,PDHPYRPYR+K_10,PYRNADNAD+K_10,NAD	Vmax,PDH = 1.87 × 10−1 K10,NAD = 10−8 K10,PYR = 1.48 × 10−2
11	V_CS=V_max,CSOAAOAA+K_11,OAAACCoAACCoA+K_11,ACCoA	Vmax,CS = 5.00 × 10−5 K11,ACCoA = 10−10−4 K11,OAA = 10−4
12	V_CITS=V_max,CITSCITCIT+K_12,CITNADNAD+K_12,NAD	Vmax,CITS = 9.33 × 10−6 K12,CIT = 1.22 × 10−5 K12,NAD = 10−8
13	V_AKGDH=V_max,AKGDHAKGAKG+K_13,AKGNADNAD+K_13,NAD	Vmax,AKGDH = 1.14 × 10−5 K13,AKG = 2.77 × 10−7 K13,NAD = 10−8
14	V_SCOAS=V_max,SCOASSCOASCOA+K_14,SCOAADPADP+K_14,ADP	Vmax,SCOAS = 2.45 × 10−6 K14,ADP = 2.00 × 10−7 K14,SCOA = 4.77 × 10−5
15	V_SDH=V_max,SDHSUCSUC+K_15,SUCNADNAD+K_15,NAD	Vmax,SDH = 4.59 × 10−3 K15,NAD = 10−8 K15,SUC = 1.27 × 10−2
16	V_FUM=V_max,FUMFUMFUM+K_16,FUM	Vmax,FUM = 5.03 × 10−6 K16,FUM = 3.44 × 10−6
17	V_MLD=V_max,MLDMALMAL+K_17,MALNADNAD+K_17,NAD	Vmax,MLD = 4.70 × 10−6 K17,MAL = 3.01 × 10−6 K17,NAD = 10−12
18	V_ME=V_max,MEMALMAL+K_18,MALNADPNADP+K_18,NADP	Vmax,ME = 8.01 × 10−7 K18,MAL = 9.01 × 10−8 K18,NADP = 10−12
19	V_PC=V_max,PCPYRPYR+K_19,PYRATPATP+K_19,ATP	Vmax,PC = 1.23 × 10−5 K19,ATP = 10−12 K19,PYR = 1.03 × 10−5
20	V_AlaTA=V_max,ALATAPYRPYR+K_20,PYRGLUGLU+K_20,GLU	Vmax,AlaTA = 10−6 K20,GLU = 10−4 K20,PYR = 10−12
21	V_GLNT=V_max,GLNTEGLNEGLN+K_21,EGLNATPATP+K_21,ATP	Vmax,GLNT = 6.83 × 10−5 K21,ATP = 6.59 × 10−6 K21,EGLN = 1.14 × 10−8
22	V_GluT=V_max,GLUTGLUGLU+K_22,GLUADPADP+K_22,ADP	Vmax,GluT = 2.90 × 10−2 K22,ADP = 3.05 × 10−6 K22,GLU = 4.43 × 10−1
23	V_GLNase=V_max,GLNaseGLNGLN+K_23,GLN	Vmax,GLNase = 1.27 × 10−3 K23,GLN = 7.16 × 10−1
24	V_GLDH=V_max,GLDHGLUGLU+K_24,GluNADNAD+K_24,NAD	Vmax,GLDH = 1.21 × 10−6 K24,GLU = 3.51 × 10−3 K24,NAD = 10−8
25	V_ASTA=V_max,ASTAOAAOAA+K_25,OAAGLUGLU+K_25,GLU	Vmax,ASTA = 10−6 K25,GLU = 10−5 K25,OAA = 10−12
26	V_AspT=V_max,ASPTEASPEASP+K_26,EASPATPATP+K_26,ATP	Vmax,ASPT = 7.01 × 10−5 K26,ATP = 10−12 K26,EASP = 6.85
27	V_AsnT=V_max,ASNTEASNEASN+K_27,EASNATPATP+K_27,ATP	Vmax,ASNT=3.23 × 10−6 K27,ATP = 10−12 K27,EASN = 4.97 × 10−8
28	V_AsnS=V_max,ASNSASPASP+K_28,ASPATPATP+K_28,ATP	Vmax,ASNS = 1.10 × 10−5 K28,ATP = 10−12 K28,ASP = 1.46 × 10−5
29	V_Asnase=V_max,ASNaseASNASN+K_29,ASN	Vmax,ASNase = 2.32 × 10−2 K29,ASN = 7.36 × 10−1
30	V_ArgT=V_max,ARGTEARGEARG+K_30,EARGATPATP+K_30,ATP	Vmax,ARGT = 3.84 × 10−3 K30,ATP = 10−12 K30,EARG = 2.42 × 102
31	V_ASL=V_max,ASLASAS+K_31,AS	Vmax,ASL = 2.00 × 10−6 K31,AS = 10−5
32	V_Arg1=V_max,ARG1ARGARG+K_32,ARG	Vmax,ARG1=4.00 × 10−6 K32,ARG = 10−4
33	V_OTC=V_max,OTCCPCP+K_33,CPORNORN+K_33,ORN	Vmax,OTC = 1.10 × 10−6 K33,CP = 10−12 K33,ORN = 10−5
34	V_iNOS=V_max,iNOSARGARG+K_34,ARGNADPNADP+K_34,NADP	Vmax,iNOS = 10−6 K34,ARG = 10−5 K34,NADP = 10−12
35	V_ASS=V_max,ASSCTRCTR+K_35,CTRASPASP+K_35,ASPATPATP+K_35,ATP	Vmax,ASS = 3.00 × 10−6 K35,ATP = 10−12 K35,ASP = 10−5 K35,CTR = 10−5
36	V_CPS=V_max,CPSNH_4NH_4+K_36,NH_4ATPATP+K_36,ATP	Vmax,CPS = 2.43 × 10−6 K36,ATP = 10−12 K36,NH4 = 4.13 × 10−8
37	V_PrAMP=V_max,PrAMPR5PR5P+K_37,R5P	Vmax,PrAMP = 6.11 × 10−7 K37,R5P = 8.90 × 10−10−7
38	V_CK=Vmax,CKfPCrPCr+K_38,PCrADPADP+K_38,ADP−Vmax,CKrCrCr+K_38,CrATPATP+K_38,ATP	Vf max,CK = 10−6 Vr max,CK = 10−7 K38,ADP = 10−12 K38,ATP = 10−12 K38,PCr = 10−5 K38,Cr = 10−7
39	V_NADPHox=V_max,NADPHoxNADPHNADPH+K_39,NADPH	Vmax,NADPHox = 1.39 × 10−4 K39,NADPH = 6.49 × 10−7
40	V_leak=V_max,leakNADHNADH+K_40,NADH	Vmax,leak = 2.95 × 10−3 K40,NADH = 10−5
41	V_Resp=V_max,RespNADHNADH+K_41,NADHADPADP+K_41,ADP	Vmax,Resp = 5.00 × 10−5 K41,ADP = 10−12 K41,NADH = 10−5
42	V_ATPase=Vmax,ATPasefATPATP+K_42,ATP−Vmax,ATPaserADPADP+K_42,ADP	Vf max,ATPase = 1.29 × 10−3 Vr max,ATPase = 8.51 × 10−1 K42,ADP = 5.11 × 10−2 K42,ATP = 2.04 × 10−5
43	V_AK=Vmax,AKfAMPAMP+K_43,AMPATPATP+K_43,ATP−Vmax,AKrADPADP+K_43,ADP	Vf max,AK = 1.55 × 10−5 Vr max,AK = 2.31 × 10−6 K43,ADP = 2.68 × 10−8 K43,AMP = 8.31 × 10−10 K43,ATP = 6.18 × 10−8
44	V_NAT=V_max,NATR5PR5P+K_44,R5PX5PX5P+K_44,X5PEGlnCr+K_44,EGlnATPATP+K_44,ATP	Vmax,NAT = 1.81 × 10−7 K44,R5P = 10−12 K44,X5P = 1.12 × 10−5 K44,EGln = 1.02 × 10−2 K44,ATP = 5.12 × 10−6
45	V_NHG=Vmax,NHGfNADNAD+K_45,NADATPADP+K_45,ATP−Vmax,NHGrNADPCr+K_45,NADPADPADP+K_45,ADP	Vf max,NHG = 4.19 × 10−7 Vr max,NHG = 1.84 × 10−8 K45,ADP = 10−6 K45,ATP = 1.05 × 10−5 K45,NAD = 10−6 K45,NADP = 10−6
46	V_growth=V_max,growthG6PG6P+K_46,G6PR5PR5P+K_46,R5PEGLNEGLN+K_46,EGLNATPATP+K_46,ATP	Vmax,growth = 4.23 × 10−3 K46,EGLN = 3.20 × 10−1 K46,G6P = 10−12 K46,R5P = 10−12 K46,ATP = 10−12

Vmax in mmol.10−6 cells.h−1 and and Km in mmol.10−6 cells and in mM for intracellular and extracellular concentrations respectively.

Table 3

Metabolites initial concentration.

Metabolite	Initial concentration		Metabolite	Initial concentration
ACCoA	1.54 × 10−8	mmol 10−6 cells	F6P	10−7	mmol 10−6 cells
ADP	3.74 × 10−7	mmol 10−6 cells	FUM	10−6	mmol 10−6 cells
AKG	6.85 × 10−7	mmol 10−6 cells	G6P	10−8	mmol 10−6 cells
ALA	10−6	mmol 10−6 cells	GAP	4 × 10−8	mmol 10−6 cells
AMP	8.25 × 10−8	mmol 10−6 cells	Gln	4 × 10−3	mmol 10−6 cells
ARG	1.4 × 10−3	mmol 10−6 cells	Glu	5.625 × 10−4	mmol 10−6 cells
AS	10−8	mmol 10−6 cells	MAL	4.1 × 10−7	mmol 10−6 cells
ASN	6.45 × 10−4	mmol 10−6 cells	NAD+	10−6	mmol 10−6 cells
ASP	3.83 × 10−4	mmol 10−6 cells	NADH	10−6	mmol 10−6 cells
ATP	2.57 × 10−6	mmol 10−6 cells	NADP+	2 × 10−6	mmol 10−6 cells
CIT	2 × 10−7	mmol 10−6 cells	NADPH	1.26 × 10−7	mmol 10−6 cells
CO2	2 × 10−5	mmol 10−6 cells	NH4	0.1287	mmol 10−6 cells
CP	10−8	mmol 10−6 cells	OAA	1.5 × 10−8	mmol 10−6 cells
Cr	3.9 × 10−5	mmol 10−6 cells	ORN	10−8	mmol 10−6 cells
CTR	10−8	mmol 10−6 cells	PCr	2 × 10−5	mmol 10−6 cells
EARG	0.669	mmol L−1	PEP	1.25 × 10−7	mmol 10−6 cells
EASN	0.313	mmol L−1	PYR	3.97 × 10−5	mmol 10−6 cells
EASP	0.184	mmol L−1	R5P	10−7	mmol 10−6 cells
EGLC	9.843	mmol L−1	SCoA	4 × 10−7	mmol 10−6 cells
EGLN	2.07	mmol L−1	SUC	3.3345 × 10−6	mmol 10−6 cells
EGLU	0.27	mmol L−1	X5 P	9.25 × 10−8	mmol 10−6 cells
ELAC	0.5033	mmol L−1	X	0.2	106 cells mL−1
ENH4	0.12	mmol L−1

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. and are respectively the experimental data and simulated values for each state variable and time , and the weight is the inverse of the variance of the experimental data for each state variable, .

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.

Table 4

Model sensitive parameter optimal values and confidence intervals.

Model parameter	Optimal value minimizing simulation error	95% Confidence intervals
		Lower bound	Upper bound
Vmax,HK	5.08 × 10−4	4.92E-04	5.79E-4
Vf max,PGI	3.53 × 10−4	1.62E-04	4.43E-04
Vmax,G6PDH	7.62 × 10−5	6.12E-05	7.38E-05
Vmax,EP	6.74 × 10−5	3.81E-05	1.09E-04
Vmax,TK	3.37 × 10−5	1.39E-05	4.28E-05
Vmax,PFK	3.42 × 10−4	3.12E-04	3.96E-04
Vmax,PK	9.79 × 10−4	9.70E-04	1.00E-03
Vmax,CITS	9.33 × 10−6	9.06E-06	1.14E-05
VAKGDH	1.14 × 10−5	1.14E-05	1.30E-05
Vmax,FUM	5.03 × 10−6	4.61E-06	5.10E-06
Vmax,PC	1.23 × 10−5	1.06E-05	2.15E-05
Vmax,ASNS	1.10 × 10−5	1.10E-05	1.10E-05
Vmax,ASNase	2.32 × 10−2	9.97E-03	2.79E-02
Vmax,OTC	1.10 × 10−6	9.77E-07	1.28E-06
Vmax,iNOS	10−6	8.90E-07	1.16E-06
Vmax,ASS	3.00 × 10−6	2.94E-06	3.07E-06
Vmax,PrAMP	6.11 × 10−7	6.02E-07	6.19E-07
Vf max,AK	1.55 × 10−5	1.53E-05	1.56E-05
Vr max,AK	2.31 × 10−6	2.27E-06	2.49E-06
Vf max,ATPase	1.29 × 10−3	8.38E-04	1.53E-03
Vmax,NADPHox	1.39 × 10−4	1.29E-04	1.35E-04
Vmax,growth	4.23 × 10−3	4.23E-03	4.42E-03
K1,ATP	1.08 × 10−5	1.05E-05	1.38E-05
K4,R5P	10−6	3.19E-07	2.00E-06
K5,R5P	10−6	4.95E-08	1.42E-06
K28,ASP	1.46 × 10−5	1.41E-05	1.51E-05
K29,ASN	7.36 × 10−1	6.91E-01	7.82E-01
K35,CTR	10−5	9.21E-06	1.07E-05
K37,R5P	8.90 × 10−7	8.66E-07	9.15E-07
K42,ATP	2.04 × 10−5	1.03E-05	2.55E-05

(Vmax in mmol.L−1.h−1 and and Km in mM).

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.

Figure 2

Simulated (line) and experimental data (dots) for BM-derived MDSCs. (A) Cell growth and nutritional profile. (B) Glycolysis. (C) TCA cycle. (D) Cell energetics. Error bars are calculated based on independent triplicates. Cell density is in 106 cells.mL−1 and intracellular and extracellular concentrations are in mmol.10−6 cells and mM, respectively.

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.

Figure 3

Sensitivity analysis of the model. Analysing most sensitive (A) parameters and (B) initial conditions. The computed value is log WSSRES and the distance is among −10%, −5%, + 5% and + 10%.

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 () and the maximal ATPase () reaction rate as well as its affinity constant () 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.

Figure 4

Time-evolution of glycolysis. (A) Glycolytic flux () in . (B) Glycolysis contribution to ATP production . (C) Comparison of pyruvate producing () and consuming () fluxes. (D) Ratio of pyruvate entering the TCA cycle (). 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.

Figure 5

Pentose phosphate pathway and glutaminolysis modulation in BM-derived MDSCs. (A) G-6-P branch point (). (B) Recirculation of glucose-derived carbon from pentose phosphate pathway into glycolysis (). (C) NADPH-to-NADP ratio. (D) Contribution of L-glutamine to TCA cycle intermediates replenishment (). All results obtained from the same model simulation than that of Fig. 2.

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 13C-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 (O2 −), 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.

Figure 6

Energetics supporting MDSCs immunosuppressive activity. (A) Percentage of NADPH consumed by L-arginine metabolism, iNOS and ARG1 activities and L-arg uptake and endogenous synthesis. (B) Percentage of ATP consumed by L-arginine metabolism. All results obtained from the same model simulation than that of Fig. 2.

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 implementing the metabolic network. This first modelling work on MDSCs thus brings a basic framework for the study of the immunosuppression phenomenon, with the most sensitive parameters as keys of further model development.

Supplementary information