GEM-Based Metabolic Profiling for Human Bone Osteosarcoma under Different Glucose and Glutamine Availability

Cancer cell metabolism is dependent on cell-intrinsic factors, such as genetics, and cell-extrinsic factors, such nutrient availability. In this context, understanding how these two aspects interact and how diet influences cellular metabolism is important for developing personalized treatment. In order to achieve this goal, genome-scale metabolic models (GEMs) are used; however, genetics and nutrient availability are rarely considered together. Here, we propose integrated metabolic profiling, a framework that allows enriching GEMs with metabolic gene expression data and information about nutrients. First, the RNA-seq is converted into Reaction Activity Score (RAS) to further scale reaction bounds. Second, nutrient availability is converted to Maximal Uptake Rate (MUR) to modify exchange reactions in a GEM. We applied our framework to the human osteosarcoma cell line (U2OS). Osteosarcoma is a common and primary malignant form of bone cancer with poor prognosis, and, as indicated in our study, a glutamine-dependent type of cancer.


Introduction
A widely recognized hallmark of cancer is reprogramming of cellular metabolism, which aims at promoting rapid cell proliferation and long-term maintenance, thus facilitating the uptake and conversion of nutrients into biomass [1][2][3][4][5][6]. Determining the metabolic requirement of proliferating cancer cells in order to modulate their metabolism might be a key factor to improve cancer treatment [7]. The genetic alterations in oncogenes and tumor suppressor genes promote cancer reprogramming cellular metabolism and supporting tumorigenesis [1,2,8]. Therefore, the energy and nutrient requirements necessary for rapid proliferation can be realized. However, genetic alteration is not the only cell-intrinsic determinant of cancer cell metabolism; the second is the origin of the cell. All cancers begin when one or more genes in a cell mutate [9]. Therefore, there is not one cancer-type cell, and each cancer cell combines the metabolic features that come from the origin tissue/organ and from the mutation(s).
Understanding the regulation of the metabolic pathways by genetic factors provides new insights into the required metabolic aspect of tumorigenesis and can also potentially provide therapeutic metabolic targets. Cancer cell metabolism similar to any other cell in living organisms is also determined by cell-extrinsic factors, including interactions with the environment [7,10]. There are numerous environmental factors that can affect cancer cell metabolism, where in particular diet, which affects nutrient availability, is one of the most important determinants [1,2,7,11,12]. Many metabolic signatures are similar across different kinds of cancer cells, such as the Warburg Effect, where cells prefer aerobic glycolysis over oxidative phosphorylation, followed by lactic acid fermentation, even in the presence of oxygen and functioning mitochondria [3,6]. However, it does not mean that all cancer cells are addicted to glucose and can not survive without it.
Here, we aim to contribute to the recently growing understanding of how certain cancer types respond to various nutrient availability and how the access to and utilization of nutrients by cancer cells affects the rate of proliferation. In this study, we showcase our approach on the human osteosarcoma cell line (U2OS); however, it could be further applied to other carcinoma and sarcoma cell lines. We provide a tool based on a genome-scale metabolic model (GEM) [13][14][15][16][17] for approaching and testing hypotheses on how nutrient availability affects cancer cell metabolism and progression. Furthermore, opposing the modulation of genetics, modulating the access to nutrients seems to be more feasible and safe. However, first the metabolic requirements of a particular cancer type should be defined. In order to determine the nutrient requirement, we propose to utilize GEMs, transcriptomic and nutrient availability data. Moreover, we investigate nutrient requirements of cancer cells in different conditions with altered availability of nutrients, and we examine its influence on transcription of metabolic genes. Our results are in line with the experimental studies carried out by other groups. It has been shown that due to the importance of glutaminolysis, most cancer cells are sensitive to glutamine deprivation and are resistant to glycolytic inhibition [18,19].
In this paper, we focused on the human osteosarcoma cell line (U2OS). Osteosarcoma is a common and primary malignant form of bone cancer. It generally affects children and adolescents, where it represents the eighth most common form of childhood cancer. With its poor prognosis, it is the second most important cause of death related to cancer in both children and adolescents [20]. This paper has a multidisciplinary character, ranging from genomics and metabolism to in silico modeling. Therefore, we state the following research questions:

1.
Do nutrient availability changes affect metabolism of the U2OS cells? 2.
Do changes in glucose and/or glutamine availability impact the growth rate and metabolic requirement of U2OS cells? 3.
Does simultaneous integration of the gene expression and nutrient availability data into a GEM provide accurate prediction of cell growth and metabolite consumption/production?
To achieve our goals we performed integrated metabolic profiling of the U2OS cell line that includes experiments and genome scale modeling. By gathering transcriptomic data from these cells in varying microenvironments, the effect of nutrient availability on cancer cell metabolism was studied. We focused on the questions whether this cell line is glucosedependent, and whether U2OS cells require glutamine in order to grow and survive. Finally, we addressed the question whether diet (using altered nutrient availability) affects the transcription level of the cells and can alter the expression of the metabolic profile. As a modeling platform we used the genome-scale metabolic models [13][14][15][16][17] that can be used to predict metabolic flux values for metabolic reactions using optimization techniques, such as Flux Balance Analysis (FBA) [16,17] and Flux Variability analysis (FVA) [21], which uses linear programming. GEMs allow to describe an entire set of stoichiometry-based, mass-balanced metabolic reactions with gene-protein-reaction (GPR) associations.
In our study, we utilized the most comprehensive human metabolic network model to date, Recon3D [22], and a novel method, named Gene Expression and Nutrients Simultaneous Integration (GENSI) [23], to create cell-specific models. GENSI translates the relative importance of gene expression and nutrient availability into the fluxes. Further, we used our models to predict the influence of glucose and glutamine availability on the biomass synthesis flux and cancer-related metabolic behaviors such as the consumption of glucose and glutamine and the production of lactate. Our integrated study that includes experiments on the cell line and the computational simulation on RECON3D-based metabolic reconstruction shows that the U2OS is independent of the glucose level, and the presence of glutamine is essential for cell viability.
The contribution of this paper is five fold:

1.
We provide an integrated GEM-based metabolic profiling framework for identification of the nutrient requirements that necessitates using RNA-seq data specific for each condition.

2.
We prove that genome-scale models like RECON3D provide a great platform for integration of transcript and nutrient availability data in order to create specific models for studying cell metabolism.

3.
We prove that translation of the relative importance of gene expression and nutrient availability data into the fluxes based on observed experimental feature(s) is a reliable method to build specific models.

4.
We identified the impact of glucose and glutamine availability on the metabolism and proliferation rate of the U2OS cells.

5.
We identified that the U2OS cell line is a glutamine-dependent cancer type.

GEM-Based Metabolism Profiling Methodology
To study the influence of nutrient availability on the growth rate and the metabolic behavior, we used the U2OS cell line and seven Dulbecco's Modified Eagle Medium (DMEM) derived media containing different concentrations of glucose and glutamine (or Lalanine-glutamine). The U2OS cell line consists of epithelial adherent cells and exhibits a fast growth rate. We were interested in the growth, nutrient requirements, and expression of the metabolic genes of the U2OS cells in optimal and not-optimal growing conditions; thus, the seven conditions were chosen (see Table 1 in Section 4).
We cultured them altering the metabolic microenvironment until the growth rate achieved a steady state. Further, we performed the proliferation assays and isolated mRNA from the cells and metabolites from the medium. We performed RNA-seq analysis using the NextSeq 550 Sequencing System (Illumina) and generated transcriptomic data after proper quality control and filtration. The RNA sequencing data were then processed (see Section 4) and used for creation of the specific GEM models. Since GEMs include known functions of protein-encoding genes, they can be used as platforms for analyzing mRNA expression data to elucidate how changes in gene expression impact cell metabolism and phenotypes [24][25][26][27][28][29][30]. Here, we additionally integrated the nutrient availability data by constraining the bounds of the exchange reactions. Our supervised approach that integrates both RNA-seq and nutrient availability data is based on the recently published GENSI [23] method that utilizes the Reaction Activity Score (RAS) approach for mapping the transcriptome into a metabolic network [25]. GENSI creates GEM-based specific models by integrating both data into the fluxes based on observed feature(s) that distinguish conditions. The proposed integrated metabolic profiling framework is schematically depicted in Figure 1. , gene expression and nutrient availability data and one feature of the cells that distinguishes cells cultured in different conditions. In the first step of the proposed approach, similarly to Gene Expression and Nutrients Simultaneous Integration (GENSI) [23], the RNA-seq is converted into Reaction Activity Score (RAS) [26], which is defined for any reaction associated with a gene in the GEM as a function of the expression of the genes encoding for the subunits and/or the isoforms of the associated enzyme(s). Nutrient Availability (NA) is converted to Maximal Uptake Rate (MUR), which is defined for an exchange reaction j in GEM and describes the rate of the maximum possible uptake over the time for substances available for cells. NA is limited by a composition of media used in our experiments. Obtained GEM-RAS-MUR models are then used for metabolic profiling using optimization methods such as FBA [16,17] and FVA [21].

Formulation of the Specific GEM Model
A central challenge in understanding and treating cancer comes from its multi-scale nature. To gain insight into the impact of the access to nutrients on metabolism and the rate of biomass synthesis, we integrated both transcriptomic and nutrient availability data into the most extensive genome-scale model of human metabolism RECON3D [23] and studied their meaning using the FBA [16,17,21] approach.
In the GENSI framework the flux bounds of the reactions that are associated with genes are set proportionally to RAS scores that are computed based on RNA-seq data. In our work, the flux bound v i of the reaction i is equal to RAS i multiplied by a factor (v i = α × RAS i ), where α is the factor independent of the reaction identity. When reaction i is reversible, its forward flux is bounded by v i and its backward flux is bounded by −v i . When the reaction i is irreversible, the flux bound in the impossible direction remains zero, and the flux bound in the possible direction is set to v i . While exchange fluxes are correlated with relative estimates of consumption and/or secretion rates. Nutrient availability data were converted into a maximum uptake rate, which describes the rate of the maximum possible uptake over the time for substances available for the model. Here, we defined MUR j for the exchange reaction j in the GEM as the absolute value of a difference in concentration of the substance s j in an extracellular environment over 48 h, for more details we refer to Section 4.
Following the GENSI framework, we found the value of the factor α by computing the steady-state flux pattern for the maximal biomass synthesis flux for various combinations of RAS and MUR within GEMs. Possible values of α were between 0.000001 and 1, see Figure 2.
For the factor value equal to 0.00042, we observed that predictions for limitations imposed by nutrient availability and gene expression, i.e., GEM-RAS-MUR, matched almost ideally maximal biomass synthesis fluxes corresponding to the proliferation rate observed in experiments (Figures 2c and 3c). Interestingly, the growth rate does not depend on the glucose concentration in the presence of 4 mM of glutamine (or L-alanine-glutamine). In the absence of glutamine, the cells were unable to survive and started dying after hours of culturing ( Figure 3e). Since cells could not proliferate without glutamine, we were not able to create GEMs for Glc10Gln0 and Glc2Gln0 due to the lack of the RNA-seq data. We then built GEM-RAS(Glc25Gln4)-MUR models using RNA-seq data from cells cultured in Glc25Gln4 medium, which was our starting condition. Computed biomass fluxes across seven different nutrient availability data (Figure 3d) clearly showed that without glutamine the U2OS cells had the least chance of survival. The biomass flux for the cases without glutamine was much smaller (Figure 3d). The results of this prediction are consistent with the viability of cells and show that the presented framework is a good approach to predict survival of cells in new conditions.
To examine whether only gene expression level has impact on the biomass synthesis, we formulated five RECON3D-based models via constraining reactions associated with genes as a function of RAS without changing exchange reactions (GEM-RAS). We did observe the impact of mapping the gene expression; however, there was no value of the factor α for which predicted biomass fluxes and growth rates were correlated (Figure 3a). This result seems to confirm that genetics is not the only determinant of the metabolic phenotype.
Next, we decided to test whether integration of only nutrient availability changes the prediction of the fluxes (GEM-MUR). We built seven GEM-MUR models via constraining exchange reactions by MUR calculated for each condition. We did not change the bounds of reactions associated with genes. As we expected, the results again did not meet the experimental observations, leading even to an experimentally unobserved situation of a high growth rate for Glc10Gln0 and Glc2Gln0 (Figure 3b). The predictions of biomass flux for GEM-RAS-MUR models were aligned with our proliferation assay (compare Figure 3c,e); the biomass flux proportional to the reciprocal doubling time was glucoseindependent and decreased as the amount of glutamine changed. for conditions Glc10Gln0 and Glc2Gln0 because cells were not able to survive and the RNA could not be isolated. ** For Glc10Gln0 and Glc2Gln0 the doubling time was not available because the growth rate was less than or equal to zero.

Flux Variability Analysis to Investigate the Nutrient Requirement
To test the prediction of the GEM-RAS-MUR models we performed FVA analysis for allowed exchange reactions while supporting biomass production rate. We further explored the range of uptake (if negative) or secretion (if positive) fluxes of glucose, glutamine and lactate by plotting the range of fluxes for corresponding exchange reactions (Figure 4).
Our main observation was that glutamine uptake was essential for all conditions. The minimum flux of the glutamine exchange reaction was always less than zero ( Figure 4c); however, for the condition Glc2Gln0.5 the value was relatively small. In medium Glc2Gln0.5 a very low amount of glutamine was required to achieve optimal growth. Whereas the minimum uptake of glucose was much less than its availability and was equal to 2 mM/48 h for most conditions (Figure 4a In silico predictions were then verified experimentally by measurement of the glucose, glutamine and lactate consumption. It turned out that all of the available glutamine was consumed by cells in all cases. While the excess glucose in the conditions Glc10Gln4 and Glc25Gln4 was not utilized by cells. Interestingly, the production of lactate was similar for Glc10Gln4 and Glc25Gln4, while cells cultured in Glc25Gln4 consumed two times more glucose than cells growing in Glc10Gln4 did.

Discussion
Cancer is one of the world's most serious health problems. Recently, cancer metabolism and nutrition have gradually become important topics in theoretical and clinical research. In contrast to traditional chemotherapy, where the goal is to kill cancer cells mainly via blocking the transcription, the nutrition-based approach assumes that cancer cell metabolism is influenced by metabolic constraints imposed by a diet. Perturbations to dietary compositions contribute to changes in plasma metabolite levels that, in turn, influence metabolite levels in the tumor microenvironment and thereby alter cancer cell metabolism and therapeutic responses. Nutrient levels in the tumor/cancer cell microenvironment can have profound impacts on cell metabolism, growth and drug sensitivity [7,[10][11][12].
Here, we follow this line of research and propose a framework for exploring the impact of the diet on cancer cell metabolism and progression using the cancer cell line and medium compositions as an experimental setup. We performed integrated metabolic profiling through enriching genome-scale modeling with metabolic gene expression and nutrient availability data. In our study, we focused on the metabolism of bone osteosarcoma, which is one of the most common primary malignant bone tumors [20]. As a GEM model we used Recon3D [22], which is a computational resource and functions as the most extensive, updated and expanded human metabolic network model. With its efficiency to connect genes to biochemical pathways, it can emphasize the promising prospects of structural analysis in genome-scale models for identifying genes related to disease and the potential for developing new treatment, biomarkers or drug repurposing. Moreover, as we recently showed, it could also be utilized as a platform for integration of the cell-intrinsic factor (gene expression level) and cell-extrinsic factor (nutrient availability) and contribute to the understanding of the impact of nutrient availability on cancer metabolism [22].
Here, we proposed a general framework of the integrated metabolic profiling that incorporates information about gene expression and nutrient availability to a GEM that could be further used to accurately predict cell growth rates and nutrient/metabolite consumption/production. For this purpose, we utilized the GENSI framework to scale reaction bounds using RAS, and then MUR to introduce information about nutrient availability to external reactions. In our experiments, we considered seven different conditions in which various amounts of glucose and glutamine (mM) were available, see Table 1 in Section 4. First, we noticed that in the cases of no glutamine (Glc10Gln0 and Glc2Gln0), all cells died within two weeks; that was a clear indication of the importance of glutamine availability for cell growth. Furthermore, the conditions Glc10Gln4 and Glc2Gln4, as well as the condition Glc25Gln4, displayed roughly similar and high growth rates (Figure 3e), indicating a more important role of glutamine compared to glucose in the growth of U2OS cells. This fact is further supported by the growth rates of the conditions Glc2Gln0.5 and Glc2Gln2, where Glc2Gln2 clearly presented a higher growth rate (almost as high as conditions with Gln4) in comparison to Glc2Gln0.5. We were able to answer our research questions (1 and 2) that the U2OS cell line is a glutamine-dependent cancer type, and glucose has very limited influence on cancer cell growth.
Next, we showed that our integrated metabolic profiling allowed us to obtain a GEM that provides highly accurate predictions. For this purpose, we first verified whether a combination of RAS and MUR (GEM-RAS-MUR) is crucial to match model predictions with observed quantities and compared them against a GEM with information about RAS (GEM-RAS) and a GEM with MUR (GEM-MUR). The results in Figure 3 seem to confirm that genetics or nutrient availability used separately are insufficient determinants of the metabolic phenotype. Only their combination (GEM-RAS-MUR) allowed us to obtain accurate predictions of the cell growth.
Lastly, we applied flux variability analysis to the GEM obtained within the proposed framework (GEM-RAS-MUR) to compare uptake/secretion fluxes of selected metabolites with measured quantities. In Figure 4 we confronted predicted fluxes by the FVA (Figure 4a,c,e) against the measured quantities (Figure 4b,d,f) of glucose, glutamine and lactate. Interestingly, the predictions given by the FVA closely reflect real consumption/production. This outcome provided further evidence that the proposed integrated metabolic profiling framework can be used to develop models capable of accurately imitating highly complex processes like cancer metabolism.
In conclusion, the results obtained in our experiments indicate that the GEM-RAS-MUR model for the U2OS cell line obtained within our framework could be used to accurately predict cell growth and metabolite consumption/production. This answers our last research question positively and opens new research perspectives for applying our framework to different types of cancer, but also other illnesses. The combination of transcriptomic data together with information about nutrients seems to be essential for developing accurate in silico models like GEMs. Lastly, our study also highlights the huge potential of GEMs as a platform to incorporate big data like RNA-seq.

Culturing U2OS Cancer Cells
U2OS (human osteosarcoma) cells were cultured in Dulbecco's Modified Eagle Medium (DMEM; GibcoBRL, Grand Island, NY, USA), supplemented with 10% Fetal Bovine Serum (FBS) and 1% antibiotics (penicillin and streptomycin) at various initial concentrations of glucose and glutamine/L-alanine-glutamine ( Table 1). Six of the compositions were made using Advanced DMEM that does not contain glucose, glutamine, phenol red or sodium pyruvate. Glc25Gln4 is commercially available (DMEM glutamaxTM). The cells were first cultured in DMEM glutamaxTM for two weeks until a sufficient number of cells was acquired. Then cells were split over seven experiments and cultured for three weeks to adopt. Cells cultured in conditions 3 and 4 started dying immediately.

Growth Rate/Cell Proliferation Assay
The total number of cells in the consequent supernatant was determined by hemocytometer counting. Mean growth rate was determined by counting cells in four nonoverlapping sets of sixteen corner squares selected at random.

Isolation and Measurement of the Metabolites from Medium
Samples from each medium were collected at 0 and 48 h, repeated three times for a total of three independent experiments for every medium composition. Metabolites were isolated using the protocol described in Sapcariu et al. [31] with minor modification. The samples were filtered using 0.22 µm syringe filters (BGB Analytik Vertrieb GmbH, Rheinfelden, Germany), freeze dried and stored until measurement at −80°C. Standard solutions for glutamine, glucose and lactate were prepared in the same way as well as control samples. High-performance liquid chromatography (HPLC) was performed using the HPLC-DAD RID LC-20AT Prominence (Shimadzu, Columbia, SC, USA) machine with a UV Diode Array Detector SPD-M30A NexeraX2 or/and a Refractive Index Detector RID 20A and an analytical ion-exclusion Rezex ROA-Organic Acid H+(8%) column (250 × 4.6 mm) with guard column (Phenomenex, Torrance, CA, USA) (5 mM H2SO4 in MilliQ water (18.2 MΩ), isocratic, 0.15 mL/min. flow rate). Injection volume was 15 µL (Autosampler: SIL-20AC, Prominence, Shimadzu), column oven temperature was 55°C (Column oven: CTO-20A, Prominence, Shimadzu) and the pressure was 29 bar.

Isolation and Obtaining of RNA Data
A total of 5,000,000 cells from each medium condition were taken for RNA analysis, also repeated three times for a total of three independent experiments for every cell condition. The RNA extraction was done using the TRI Reagent Protocol from Sigma Aldrich (https://www.sigmaaldrich.com/technical-documents/protocols/biology/trireagent.html), after which all samples were treated with DNase using the TURBO DNAfree Kit from Invitrogen. Nanodrop was used for determining RNA concentrations when needed. With the Agilent 2100 Bioanalyzer System the quality of the RNA was checked, as well as the RNA Integrity Number (RIN), using RNA Pico Chips. rRNA depletion and library preparation for sequencing was completed with the NEBNext protocol ('Protocol

RNA-Seq Data Processing
RNA-seq data were available in .fastq file format. There were three measurements per condition. The RNA-seq data were processed according to the following pipeline: 1.
Sequences were mapped against the human reference genome using Burrows-Wheeler Aligner (BWA; version 0.7.17, https://sourceforge.net/projects/bio-bwa/). The outputs of this step were files in .sam format.

Mapping Transcriptomic Data to the Fluxes with Recon3D
Mapping transcriptomics data to the fluxes with Recon3D included two steps: (1) conversion of RNA-seq into RAS, (2) RAS to fluxes of reactions with gene-protein-reaction association rules.
RNA-seq data were converted to RAS, for each condition and each reaction, based on protocol outlined by Graudenzi, A. et al. [26] according to gene-protein-reaction association rules (GPRs). The formula includes AND and OR logical operators. AND rules are employed when distinct genes encode different subunits of the same enzyme, i.e., all the subunits are necessary for the reaction to occur. The AND operator is implemented as the minimum operator [23,25,26]. OR rules describe the scenario in which distinct genes encode isoforms of the same enzyme, i.e., either isoform is sufficient to catalyze the reaction. The OR operator is implemented as the summation operator [23,25,26].
In the next step, RAS i is converted to a lower bound v i of the reaction i having GPR association rules as follows: where v is the vector of fluxes through all reactions having a gene rule, v il and v iu are the vectors of lower and upper bounds on the reaction i.

Mapping Nutrient Availability (NA) to Fluxes with Recon3D
Mapping NA to fluxes with Recon3D includes two steps: (1) conversion of NA to Maximum Uptake Rate (MUR), (2) conversion of MUR to fluxes of exchange reactions. We define a MUR j for each exchange reaction j based on NA data for each condition. In the next step MUR j is converted to a lower bound v jl of exchange reaction i as follows: where v i is the vector of fluxes through exchange reactions with the environment of the system, v jl and v ju are the vectors of lower and upper bounds on these fluxes, and MUR is a Maximum Uptake Rate that is defined as the rate of the maximum possible uptake over the time for each substance available to the model. In our case, MUR is defined for the exchange reaction j in a GEM as the absolute value of the difference in the concentration of the substrate S j in extracellular environment (|[S jF ] − [S jI ]|) over the 48 h.

Simulations
In the simulations of maximal biomass flux we applied the computational technique of Flux Balance Analysis (FBA) [16] to the human genome-wide metabolic map Recon3D [22] using the COBRA toolbox in MATLAB and Python [32,33]. FBA solves the following linear programming problem: where S is the stoichiometry matrix indicating how many molecules of each metabolite are produced or consumed in each reaction, v is the vector of fluxes through all reactions including exchange reactions with the environment of the system, v l and v u are the vectors of lower and upper bounds on these fluxes, respectively, and c is a vector of weights generating the linear combination of fluxes that constitutes the objective function Z.
In the simulations of maximal and minimal fluxes for defined exchange reactions we applied the computational technique of Flux Variability Analysis (FVA) [21]. FVA is used to find the minimum and maximum flux for reactions in the network while maintaining some state of the network, e.g., supporting maximal possible biomass production rate. FVA entails the following linear programming problem: where Z 0 = w v 0 is the optimal solution to the FBA problem with biomass reaction as the objective function, w represents the biomass objective (vector of weights generating the linear combination of fluxes that constitutes the objective function Z), v is the vector of fluxes through all reactions, and γ is a parameter that controls whether the analysis is done with respect to suboptimal network states (0 ≤ γ ≤ 1) or to the optimal state (γ = 1).

Data and Source Code Availability
All related data sources such as gene featureCounts, RAS scores, and GEM models can be found at https://github.com/e-weglarz-tomczak/GEM-RAS-MUR. Data Availability Statement: The data for this paper is available at: https://github.com/e-weglarztomczak/GEM-RAS-MUR.