Nature and Nurture integrated through an adjustable Flux Balance Analysis: The impact of cell-nutrition changes on the Warburg effect in hepatocellular carcinoma

Cancer cells are often home to profound metabolic specialization, such as the Warburg effect wherein cancer cells exhibit high glucose uptake with conversion to lactate in the presence of oxygen. To understand the mechanisms underlying the Warburg effect in hepatocellular carcinoma and to investigate whether and how both the essential nutrient concentrations and the cells’ gene expression, impact on this metabolic alteration, we integrated computational genome-wide metabolic modeling with an experimental study on cell lines. We converted transcriptome sequencing (RNA-seq) data from two hepatocellular carcinoma lines HuH7 and PLC/PRF/5 to flux bounds in the most comprehensive human metabolic network model to date, Recon3D. A new method was developed to enable changes in medium concentrations to impact Flux Balance Analyses (FBA) of genome wide metabolic maps. By integrating the two methods, the cells’ metabolic behavior could be studied due to applying a new Nature Nurture Scaling (NNS) factor, which assesses the relative importance of gene expression (‘Nature’) and medium composition (‘Nurture’). To examine whether the medium concentrations of glucose and glutamine affect the rate of aerobic glycolysis and growth rate in vitro, we further cultured the same two hepatocellular carcinoma lines on media that contained different amounts of glucose and glutamine. Our integrated approach where expression level (nature) and the environment (nurture) work together, led us to conclude that proliferation, glucose consumption, and lactate production are associated with the presence of glucose, but do not necessarily increase with its concentration when the latter exceeds the physiological concentration. There was no such association with the presence of glutamine. The observed dependencies on glucose concentration could hereby be understood in terms of a balance between Nature and Nurture, as could be the lack of the cells’ response to glutamine. Author summary Metabolism is directly involved in many human diseases including cancer, and indirectly in virtually all, because disease causes metabolic changes that can accompany or even affect etiologies and be read as biomarkers. The best-known metabolic abnormality in cancer cells is an increased glycolysis followed by lactic acid production even in the presence of oxygen and fully functioning mitochondria, a process known as the Warburg Effect. So-called genome-scale metabolic reconstructions integrate all known metabolic reactions occurring in an organism into a single map. Together with a mathematical method for simulating the optimal balance of metabolic fluxes, this holds promise for studying the mechanisms by which networks control life. We combined this genome-scale metabolic modeling approach Węglarz-Tomczak et al. Understanding Warburg fluxes, page: 3 with in vitro experiments to investigate whether the behavior of cancer cells is determined by their nutrition or/and the expression of their genes. We hereto developed a new method that integrates the genome (nature) and the environment (nurture) and identified the influence of cell-nutrition changes on the Warburg effect in hepatocellular carcinoma. Highlights: 1. Cell nutrition and gene expression can now be dealt with simultaneously by FBA through a newly developed Nature Nurture Scaling (NNS) factor. 2. The physiological concentration of glucose is sufficient for the production of lactate by Huh7 and PLC. 3. Huh7 and PLC growth rates are independent of glutamine in media containing glucose. 4. On glutamine alone Huh7 and PLC growth rates are much lower than on glucose alone.


Introduction
Hepatocellular carcinoma (HCC) is among the most common primary liver malignancies [1] and a frequent cause of cancer death [2]. It represents the fifth most common cancer worldwide and over 500 000 new cases are diagnosed per year [3,4]. Hepatotropic viruses such as HBV, HCV, and hepatitis D virus (HDV), and excessive alcohol intake leading to chronic liver disease and cirrhosis, are the most common cause of HCC worldwide [5,6].
There are other relevant risk factors such as tobacco, obesity, diabetes, metabolic syndrome, and selected aspects of diet [5]. Several factors contribute to the gloom of HCC prognosis: 1) late diagnosis due to lack of symptoms at early stages, 2) recurrence and resistance to chemotherapy, 3) large heterogeneity, and 4) HCCs perfect environment of healthy liver cells providing nutrients to their surroundings. Notwithstanding therapies such as liver resection, transplantation and ablation, the prognosis of HCC is still very poor [7] with five-year survival rates in the 20% range in the United States [8].
Altered metabolism in cancer has been linked to increased tumorigenicity and resistance to chemotherapeutics. The metabolic signatures of cancer cells result from biochemical/metabolic reprogramming (switching) that promote fast growth, survival, proliferation, and long-term maintenance, and are thereby selected for [9]. In most cancer Węglarz-Tomczak et al. Understanding Warburg fluxes,page: 4 cells the rate of glucose uptake has increased and part of glycolysis-derived pyruvate is diverted to lactate [10][11][12][13], which produces much less ATP per glucose than its oxidation to carbon dioxide would. This metabolic signature, known as the Warburg effect [12], enables dividing cells to satisfy anabolic needs for biomass production and is accompanied by a suppression of apoptotic signaling [14,15,16]. The high glucose consumption during the Warburg effect also provides a higher production of reduced nicotinamide adenine dinucleotide phosphate (NADPH 2 ) through the pentose phosphate pathway, which provides electrons to cell proliferation [2,17]. Cancer cells tend to persist in converting much glucose to lactate [10][11][12][13], while most normal cells import the glycolysis-derived pyruvate into the mitochondrial matrix. Here the pyruvate is converted to acetyl-coenzyme A, which is oxidized by mainly NAD+ in the Krebs cycle to CO 2 [10,11]. The Warburg effect [12] correlates with tumor aggressiveness and poor patient prognosis in many tumor types, including HCC [14,[18][19][20][21]. Many cancer cells also have an increased uptake of glutamine [22][23][24][25][26]. The partial catabolism of this glutamine to lactate by cancer cells has been called the WarburQ effect [27]. Some rapidly proliferating cells are particularly dependent on glutamine, and undergo necrosis upon glutamine depletion [25], possibly because glutamine is their main source of nitrogen essential for the synthesis of the amino acids and nucleotides required for cell proliferation.
The excess secretion of lactate would appear to be an inefficient use of cellular resources: each excreted lactate molecule wastes three carbons that might otherwise be utilized for either ATP production [12]. However, the human organism is a perfect host for cancer cells. Lactate produced by cancer cells is picked up from the bloodstream and reconverted to glucose by the liver via the Cori cycle [28], the glucose returning to the blood that then supplies cells, including tumor cells, with this glucose. Similar to lactate, ammonia produced during high glutamine catabolism can also be removed by the liver via the putative glutamine-ammonia cycle [29,30] and returned to the tumor cells as glutamine. According to "Lactate Shuttle" concepts, lactate plays a role in delivery of oxidative and gluconeogenic substrates as well as in cell signaling [31][32][33]. It includes intracellular (e.g. between the cytosolic and mitochondrial compartments) and extracellular (e.g. between working skeletal muscle and heart) lactate exchanges [31]. The difference between processes e.g. in the healthy heart and skeletal muscles of the host, and cancer cells, is lack of regulation and asocial behavior of cancer cells [27]. The host's physiology bears the burden of providing a relatively unlimited glucose supply to the cancer cells and a sink for the lactate they secrete [10a,10b]. Thereby the host delivers resources to the proliferating cancer cells at the cost of delivery of these same resources to the other organs, and increasing whole-body Gibbs energy expenditure and multiple organ failure as is highlighted in cancer cachexia. In the case of HCC, resourcing processes are simplified and the liver becomes the host in the host.
Refurbished glucose is readily used by local cancer cells and lactate produced by cancer cells is used by liver cells in gluconeogenesis without transportation through the body. The switch is not always quite irreversible: switching of glycolysis to gluconeogenesis in hepatocarcinoma cells lacking gluconeogenesis was possible after application of the corticosteroid medication dexamethasone [34].
A better understanding of the mechanistic links between cell metabolism and survival control may aid in the development of strategies towards a more specific control of HCC progression as well as other types of cancer. The genome-scale metabolic models (GEMs), that are at the core of some bottom-up systems biology approaches, constitute promising tools to study such complex diseases [35,36,37]. Using metabolic reconstructions, researchers can store and continually update information about chemical reactions in a standardized genome wide representation. Various human metabolic maps that represent the entire network of metabolic reactions in a human is known, have been published and their content has been expanded e.g. from 3311 reactions in RECON1 [35a] to 13,543 reactions in RECON3D [37]. Recon3D is not only the most comprehensive human metabolic network model to date, it also includes three-dimensional (3D) metabolite and protein structure data and enables integrated analyses of metabolic functions in humans. GEMs have already been used to predict biomarkers of inborn errors of metabolism [36,38], to identify drug targets [39], to characterize cancer metabolism [40] and to improve the understanding of microbial interactions with the host organism [41,42]. This has been achieved by a computational approach called Flux Balance Analysis (FBA), which enables the computation of optimal steady state flux distributions over genome-scale metabolic maps, optimal in terms of objective functions such as maximum growth yield. Existing computational methodology is able to calculate the implications of quantitative gene expression information for those optimal flux patterns and yields. In reality, the performance of cells may not only be a function of gene expression (which we shall here call 'Nature') but also of the metabolic environment of the cells ('Nurture'), but there is no corresponding FBA methodology yet to take nutrient concentrations into account quantitatively in such calculations.
In the present study, we developed the latter methodology and then use the GEM to predict the influence of the level and ratio of two essential nutrients (glucose and glutamine) in medium on the growth rate and metabolic behavior of cancer cells (HCC), whilst taking into account the transcriptome of those cells. We took advantage of the novel metabolic reconstruction RECON3D [37] and created RECON3D-based metabolic maps that integrate, for the first time, transcriptomic (nature) and nutrient concentration (nurture) data through a new Flux Balance Analysis approach, by using a new 'Nature Nurture Scaling factor' ('NNS'). Using this integrated approach, in silico modeling and in vitro experiments on the two hepatocellular carcinoma cell lines HuH7 and PLC/PRF/5, we identify how limitations in expression levels of the metabolic network and nutrient limitations may both cause glucose and glutamine (in)dependencies.

Cell lines
For our study we chose two HCC cell lines, i.e. the hepatitis infection-negative Huh7 and the hepatitis infection-positive PLC/PRF/5. Huh7 is a well differentiated hepatocyte-derived cellular carcinoma cell line originating in a liver tumor in a 57-year-old Japanese male. Huh7 is an immortal cell line of epithelial-like, tumorigenic cells. Its cells adhere to the surface of flasks or plates and typically grow as 2D monolayers. Also known as the Alexander hepatoma cell line, PLC/PRF/5 is a human hepatoma originally taken from the liver of a patient with primary liver cell carcinoma who was persistently infected with hepatitis B virus.
Also, these cells adhere to the surface of flasks or plates and grow as 2D monolayers [43,44]. Human HCC cells (Huh-7D12 (ECACC 01042712) and PLC/PRF/5 (ECACC 85061113) were purchased from the European Collection of Authenticated Cell Cultures (Salisbury, UK).

Microarray data
We obtained microarray transcriptomics data from the MERAV database [45] (http://merav.wi.mit.edu/) for both the Huh7 and PLC/PRF/5 cell lines. There was data available from two experiments for Huh7 and one for PLC/PRF/5. All MERAV microarray datasets were renormalized together [45]. For the two Huh7 experimental datasets we averaged the signal per gene between the two experiments.

Flux balance analysis
In all simulations in this work we apply the computational technique of Flux Balance Analysis (FBA) [47] to the human genome-wide metabolic map Recon3D [37] using the COBRA toolbox in MATLAB and Python [48,49]. FBA entails 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 considered, ɑ and β are the vectors of lower and upper bounds on these fluxes, and c is a vector of weights generating the linear combination of fluxes that constitutes the objective function Z. A flux distribution resulting from FBA therefore satisfies the requirements that each metabolite is produced at the same rate as it is consumed, that the flux boundaries are not exceeded and that the flux distribution maximizes (or minimizes) the objective function Z.

Computational media representation
Ideally, to reproduce experimental cellular conditions accurately in terms of the flux distribution, one should experimentally measure the transport fluxes of all medium components that are taken up or secreted. Using FBA one can then compute the internal flux pattern that should be optimal with respect to a certain objective. However, this is impractical for most mammalian cell line studies as there are quite a few such medium components. Therefore, we here made several assumptions that may be approximately correct under very strict conditions. We assume (i) that the concentrations of the medium components are essentially constant in time, (ii) that the uptake fluxes depend by first-order (i.e. proportional) kinetics on the known substrate concentrations in the medium, thereby producing an uptake rate that can only be reduced by changes in intracellular metabolite concentrations, and (iii) that the proportionality constant is the same for all medium components. Assumption (ii) biochemically implies that in order for the uptake rate to be reduced, something in the downstream network needs to be limiting (i.e. an enzyme concentration). This ultimately has to translate back to increased concentration of metabolites downstream of the uptake step. This is only straightforward if the concentrations are well below the K m of the transporters, but may be replaced by the options that the uptake reactions run at a V max with respect to their extracellular substrate concentration but that their substrate induces them such that the transporter expression level is proportional (again with the same proportionality constant for all substrates) to (or increasing monotonically with) the extracellular concentration, and that either this induction or the rate is only negatively influenced by intracellular concentrations. Assumption (iii) is unrealistic biochemically, but may be seen as a way to reflect the metabolic potential offered by the media. In accordance with these assumptions, we set the maximal uptake fluxes, for all medium components listed in Table 2, equal to their explicit concentrations in the medium.
In FBA there is a distinction between exchange of a metabolite and transport of that metabolite. Exchange refers to the net consumption or production of a metabolite and occurs between the extracellular compartment and the outside world. It measures how much substrate that is added to the experimental medium is taken up and used by the cell and how much product leaves the system. In contrast, transport refers only to the transfer of a metabolite between the cytosol and extracellular compartment. Technically, this distinction exists to allow some metabolites to circumvent the steady state requirement in FBA [47]: net catabolism of glucose is consistent with steady state (constant concentrations of metabolites) by allowing exchange with an external compartment to which the steady state requirement is not applied. By implication, if the concentration of a metabolite in the medium supplied to the cells is zero, then the bound (=upper limit) on inward exchange flux is zero.
This does not mean that the transport reaction of that metabolite is blocked (has bounds of zero) but it does mean that that the transport flux will be zero due to lack of substrate. If that extracellular concentration were to go up (by increasing the exchange reaction bound), the metabolite could be imported and transport flux would be possible.
We did not account for metabolic components spilling over from the added serum.
We blocked all other uptake of metabolites so that only those compounds listed in Table 2 were allowed to be exchanged in. We did not alter the Recon3D default choices of which metabolites may be net produced by the in-silico cell. This left 1559 metabolites which are allowed to be net-produced by the cell (Supplementary Excel Table 1). Recon3D also contains various so-called sink and demand reactions which serve as sources and sinks for certain metabolites allowing them to bypass the regular mass-balance. We blocked all such reactions. All uptake reactions and their maximal uptake bounds were taken the same across the six conditions with the exceptions of glutamine and glucose uptake which were varied as specified in Table 1  concentrations (colored in green) and were exempt of ammonia. All concentrations were effected in-silico as maximal uptake rates which shape and constrain the flux cone of the solutions in FBA (see text). For carbon dioxide and oxygen, a not limiting uptake bound of 1000 was taken. Oxygen was considered non-limiting because the concentration in the medium was in the order of 0.2-0.4 mM (corresponding to air saturated saline) [50] whereas the K m of cytochrome oxidase for oxygen is some 0.01 mM [51]. Carbon dioxide was considered non-limiting due to its continuous replenishment in the medium.

2.3.3.
Adding an acylgroup (Rtotal) synthesis reaction to Recon3D [37] We observed that with the default Recon3D model and given the medium definition discussed above, biomass synthesis was not possible. We traced the problem back to the triglyceride synthesis pathway where in the default Recon3D version a 'source' reaction for triglycerides is present (which allows influx into the cell independent of the presence of triglycerides in the medium; it may be noted that in Recon3D such a source reaction is called a sink reaction, by virtue of sign notation; uptake fluxes are positive, effluxes are positive [47]). Indeed, when temporarily reactivating the triglyceride source reaction (or equivalently activating the triglyceride exchange reaction and adding triglyceride to the medium) a biomass synthesis became possible in all conditions.
Biochemically, triglycerides are synthesized starting from glycerol-3-phosphate and various lipid tails esterified to CoA. Each of these lipid tails can assume any of the three positions in the triglyceride molecule. We observed that without adding triglyceride uptake (or a ditto source) to the metabolic map, the first intermediate in the pathway (lysophosphatidic acid; the monoglyceride with a phosphate on the 3 position) could not be net-produced (see the network diagram in Figure S1 and S2). Inspection of the network revealed that this was due to the need for net input of 'Rtotal', 'Rtotal2' and 'Rtotal3' groups in this pathway for which there exists no synthesis reaction in Recon3D. This problem affects biomass flux because multiple metabolites downstream of the triglyceride synthesis pathway (e.g. Phosphatidylcholine, Phosphatidylserine, and Phosphatidylethanolamine, see Table   S2), or in branches of it, are explicit components of the biomass used in Recon3D. Recon3D does have the potential to make Stearoyl-CoA, Palmityl-CoA, Oleoyl-CoA and Octadecadienoyl-CoA, but lacks the reactions to associate these with the glycerol moiety: it instead associates Rtotal, Rtotal2 and Rtotal3 to the glycerol moiety, the numbers referring to the position they take in the resulting triglyceride molecule. Recon3D worked around the ensuing problem of lack of biomass synthesis by adding a source for triglycerides. We removed such dei ex machina by forbidding source reactions and thereby came across this problem. We solved it by equating Rtotal, Rtotal2 and Rtotal3 species in Recon3D to a single pool Rtotal and by adding a pooling reaction for lipid tails: Rtotal2CoA and Rtotal3CoA respectively, could be could synthesized but these metabolites were not connected to anything else in the network. Our grouping of the synthesis of all Rtotal into a single reaction, has the advantage that a ratio may be imposed when this is known experimentally. We here took this ratio to be 1:1:1:1 This workaround allowed Recon3D to sustain a positive biomass flux on the medium discussed above. An alternative in which all Rtotal-CoA could be synthesized from any of the four above CoA esters was not here entertained. The resulting 'patched' version of Recon3D is available in the supplementary code and model archive and a list of its reactions and metabolites can be seen in Supplementary Excel Table 2.
All simulations discussed in this section may be reproduced using the scripts in the Supplementary Code Repository.

Mapping transcriptomics data to Reaction Activity Scores with Recon3D [37]
We converted the entrez gene identifiers in Recon3D to their gene symbols (e.g. 8639.1 was converted to AOC3), using the mygene module in Python (https://pypi.org/project/mygene/), in order to match them to genes represented in the transcriptomics datasets. Out of the 2248 genes in Recon3D, we were unable to match 103 to the dataset on the basis of their gene symbol alone. We then additionally searched the dataset for known gene aliases and this led us to identifying an additional 87 transcripts. 16 genes of Recon3D (see Table S1) we could not identify. These included the three mitochondrial genes encoding the 3 subunits of cytochrome c oxidase. We artificially assigned all these 16 genes a TPM score equal to the maximum TPM score of the genes we could identify for each cell line. Since we do not know whether or not these genes are expressed, we did not want to artificially block them in our model analyses. Setting them to the maximal observed TPM value guarantees that they will not be limiting in any of our analyses.
The existing RNA-sequencing methodology suffers from so-called zero-inflation [53], i.e. the lack of transcripts for genes that are in fact expressed. For data integration this is problematic since a single zero may block an entire pathway. For our dataset we did indeed observe this problem. For example, in the RNA-seq data the serinepalmitoyltransferaselong-chain-base-subunit-3 gene (SPTLC3) which catalyzes the reaction synthesizing 3dehydrosphinganine (SERPT), has a TPM of zero, and blocking this reaction (after imposing our changes to Recon3D as discussed above) singlehandedly prevents biomass synthesis.
To bypass this problem, we used a microarray dataset and calculated the ratio R of the microarray intensity divided by the genome-wide median for each gene that came associated with a TPM score of zero in the RNA-seq dataset. Then we updated such genes' TPM scores and set them equal to the calculated ratio R for the microarray dataset multiplied by the genome-wide median in the RNA-seq dataset. Below we will therefore refer to this set of zero-adjusted TPM scores as TPM*. Because the microarray dataset had scores for all genes, this removed all zeros from the dataset. We here neglected any transcriptome difference between the cell lines and experimental conditions used for the microarray and RNA-seq experiments and we assumed that the microarray data were quantitative also at low gene expression.
Recon3D contains an annotated gene-reaction coupling rule for each reaction. Using AND and OR logic this rule specifies which genes encode proteins that may help catalyze that reaction. The AND logic may be used to indicate proteins that consist of more than one subunit, or protein complexes that catalyze reactions whereas the OR logic may be used to indicate isoenzymes or alternative configurations of the protein complexes. When integrating the transcriptomics data into the map we turned these Boolean gene-reaction coupling rules into quantitative rules. Here we were inspired by the Metabolic Reaction Enrichment Analysis (MaREA) methodology [54] and the E-flux approach [55]. MaREA had been developed to compare reaction activities between patients rather than between cell lines by assigning a quantitative so-called Reaction Activity Score (RAS) to each reaction based on gene expression levels for all proteins that might be involved in the catalysis of that reaction also turning the above Boolean rules into quantitative activities. In order to do this, one needs to know the levels of the corresponding proteins and protein subunits in the cell of interest.
We used the mRNA levels as a first approximation to the corresponding protein levels, assuming that the two were proportional. We assigned to each reaction a RAS by summing over isoenzymes (for OR logic) and taking minima of subunits of a complex (for AND logic) the TPM scores for the genes coupled to each reaction. In this way, isoenzymes are thought to contribute additively to the activity of a reaction whereas lack of even one subunit of an enzyme complex can linearly bring down a reaction's activity [54,55].

Constraining fluxes based on RAS scores
To evaluate the implications of the RAS scores for all reactions in the map for the metabolic fluxes through the map we mimicked the approach used in [54,40]  Proliferation assays were conducted in 25 cm 2 T flasks, starting with a cell density of 8. genes. We extracted the RNA-seq records for 2232 out of the 2248 metabolic genes that surface both in this data set and in Recon3D (see Methods). Comparing the mRNA levels (in terms of TPM scores) for the subset of metabolic genes to the genome-wide mRNA levels, we observed that most metabolic genes exhibited a higher than average expression, with a median expression level of ~15 compared to ~1 genome-wide and a mean expression level of ~81 compared to ~39 genome-wide.
We were confronted with a well-known issue with RNA-seq data integration in metabolic models, i.e. zero expression levels, including zeros in mRNAs encoding enzymes catalyzing essential reactions or combinations of enzymes that are essential. Of the 2232 metabolic genes 262 and 310 had TPM scores equal to zero in Huh7 and PLC respectively.
Such zeros could be due to true absence or correspond to technical zeros where the gene had in fact been transcribed but was somehow not measured. Technical reasons may include inefficient cDNA synthesis due to tertiary structure formation, amplification bias, or low sequencing depth. Additionally, zeros may occur due to transcription bursting between somehow synchronized individual cells [56] or to small time-windows of expression. Given that an independently obtained set of microarray data might not suffer from quite the same problems, we assigned to genes with zero TPM scores in the RNA seq analysis, alternative TPM scores that reflected the microarray datasets (see Methods).  Figure   1A).
We first looked at the genes that did correlate: In the RNA-seq dataset (prior to correcting for zeros as discussed above), 186 metabolic genes (~8% of the 2232 genes) did not come with any transcript in either cell line. 215 metabolic genes (~10%) exhibited nonzero TPM values below the genome-wide median (~1 TPM), again in both cell lines.
Together, these constituted a common set of 401 metabolic genes that were expressed at low level. On the high-abundance side, the two cell lines shared 1422 genes (~63%) that were more highly expressed than the median gene. A 623-genes subset of these 1422 (~28%) was commonly expressed above the genome-wide mean (~39 TPM). These genes behaved in line with our expectation of metabolic similarity between these cell lines. They came however with a possibly important exception of some 425 genes that were off-diagonal in Figure 1A. 415 metabolic genes (~18%) exhibited a TPM score above the genome-wide median in both cell lines and a differential expression ratio of at least 3 (or below ⅓). The data suggests that our expectation was not quite right: there was an appreciable metabolic difference between the two cell lines. In accordance with the above, the expression ratio between the two cell lines was still mostly distributed narrowly around 1. The corresponding probability distribution was largely log-normal, but not quite: it had long 'tails' on either side, suggesting that a disproportionate number of metabolic genes were much more expressed in Huh7 than in PLC and vice versa ( Figure 1C).  Recon3D dealt with this issue qualitatively through its gene-reaction coupling rules. We used a quantitative version of these rules to assign a 'Reaction Activity Score' (RAS). The RAS for a metabolic reaction reflects the expression levels of isoenzymes and components of multi-component complexes that may catalyze that reaction (Methods and [40,[54][55]).
The MaREA data integration framework [54] put forth this approach for integrating transcriptomics data into genome-scale metabolic models based on earlier work on the Eflux approach [55]. We then asked whether the metabolic differences we found between the two cell lines would disappear when correcting for these isoenzyme and enzyme subunit issues by assigning reaction activities. In Figure 1B we correlate the RASs between the two cell lines, and panel 1D shows the distribution over the metabolic genes, of the RAS ratios between the two cell lines. The RAS correlation between the two cell lines is only a little stricter than that of the individual mRNAs. The standard deviation in the RAS ratios is ~14 compared to a standard deviation of ~177 for the TPM* ratios. The fraction of outliers outside the lognormal distribution of the ratio remains substantial however. Supplementary Excel Table 3 lists the 92 outlier reactions with a log 10 RAS ratio > 2 or < -2. RASs into metabolic map flux bounds using a method inspired by the Metabolic Reaction Enrichment Analysis (MaREA) data integration framework [54]. Even though it has been suggested that RAS scores of individual reactions are better for highlighting metabolic differences among groups of (patient) samples than flux distributions predicted based on those RAS scores, here the interest is in flux distribution predictions. The authors of [54] already applied the E-flux methodology to incorporate the RAS scores into flux bounds. We approach in that it may be inaccurate due to non-proportional translation (protein synthesis) [58], absence of corresponding steady state conditions, or due to the differences in catalytic efficiency (k cat ) between enzymes. Notwithstanding these shortcomings, the approach may provide some insight into the potential network effects of changes in gene expression levels.

A new Flux Balance
Our variant of the E-flux approach allows us to scale the RNA-seq data, in reference to a specific set of medium uptake bounds (see above).
It is a-priori unclear how large the NNS factor should be. There are two possible sources of limitation for the model: the medium composition (as represented by maximal exchange rates) and the enzyme expression levels (consequent to transcriptomics). We here wish to examine the case where the enzyme expression levels begin to impose limitations on the model output. We thereto start with high α values, where the enzyme expression levels are not limiting, and maximal biomass flux is determined by substrate concentrations in the medium. Then we decrease and hence all resulting metabolic flux bounds uniformly until we see differential effects on the predicted maximal growth rates across media conditions. Where this occurs significantly, we fix α, which thereby becomes a fitted parameter.

Accounting for medium metabolite concentrations in flux balance analysis
To predict the effects of nutritional differences in terms of all components in the medium including the various concentrations of glucose and glutamine, between our six different medium conditions (see Table 2), on the global metabolic behavior of the two cell lines we modified existing genome-scale flux balance analysis: We equated the upper bound of each uptake reaction to the corresponding metabolite concentration in the medium, setting the bound to zero if the metabolite was absent from the medium. In [54,55] (Figure 2). For either cell line Figure 2 shows that for a NNS factor in excess of 0.3, the predicted growth rates differed between media conditions, the ones at 25 mM glucose being about double those at 5.6 mM glucose. Also, a dependence on glutamine concentrations was predicted, but this dependence was smaller.
When decreasing the NNS factor to below 0.1, this medium dependence disappeared for four out of six modeled medium conditions: the same biomass flux was then predicted which was still significantly higher than that for the remaining two medium conditions at zero glucose (5 and 6). Huh7 predictions for media 1-4 converged already at higher values for than did the predictions for PLC. Huh7 was predicted to have equal or lower growth rates than PLC across all conditions. Medium 1, followed by medium 2, were predicted to yield the highest biomass fluxes for the high values of , corresponding (see above) to absence of gene-expression limitations. This reflects the model's sensitivity to carbon input for high values, since medium 1 and 2 contained the highest levels of glucose. At = 0.1, the simulations for media 1-4 yielded equal biomass fluxes which were still larger than the predicted fluxes for medium 5 and 6 which lack glucose. This indicates that by reducing the NNS factor, the model can be made more (high , hence no limitation by low transcription of metabolic genes and thereby limitation by uptake) or less (low , hence strong limitation by low transcription of metabolic genes) sensitive to variation in concentration of the growth substrate. The fact that media 1-4 converge to similar biomass synthesis flux optima for low levels of may reflect a shared limiting reaction downstream of (and at a flux bound smaller to the bound of) the exchange reaction the flux bound of which keeps monotonically decreasing with decreasing . In Figure 3 we summarize the predicted maximal biomass fluxes for the NNS factor = 0.1, which is at the transition between limitation by extracellular substrate levels and intracellular expression levels. It shows that reduction of glucose concentration does not decrease maximal biomass flux as we observed in in vitro experiments (compare Figure 3 and 6A) .
We observe that nontrivial predictions for limitations imposed by medium composition and gene expression can be computed by our new method, such as that (i) both in the absence and in the presence of glutamine the growth rate should be independent of glucose concentrations between 5.6 and 25 mM, yet decrease appreciably in the absence of glucose, (ii) the specific growth rate of Huh7 cells is lower than that of PLC cells, (iii) in the absence of glucose, the cells should be able to grow on glutamine, but (iv) growth rate on glutamine alone should be much lower than on glucose alone.  Table 1 Figure 2).

Metabolic flux potential as predicted by Flux Variability Analysis
Fluxes through intracellular biochemical reactions are not only determined by the concentrations of the enzymes catalyzing them, but also by metabolic regulation [57], i.e. by the concentrations of intracellular metabolites. FBA is oblivious of metabolic regulation other than that it philosophizes about what flux should be optimal for the cell in view of some objective. The transcriptome and extracellular-concentrations informed flux bounds that we here implemented, merely define ranges of the fluxes rather than that they precisely predict the fluxes. For the precise predictions of fluxes one needs fully dynamic models [58], but the kinetic information required for this approach is missing for mammalian cells. We hereby can only predict the ranges of fluxes that are consistent with transcriptome and extracellular metabolite concentrations and this is done by flux variability analysis.
For the NNS factor value = 0.1, indicated by the black vertical line in Figure 2, we analyzed each of the 12 models in terms of the possible ranges its production fluxes of lactate, pyruvate, ATP and CO 2 and its uptake fluxes of glucose, glutamine, oxygen and phenylalanine. Here we maintained the maximal biomass flux for each specific medium and transcriptome (i.e. the biomass fluxes listed in Figure 3, which differ between the 12 models) ( Figure 4). The results for lactate and pyruvate production are non-negative since these compounds are not in the growth medium and thus cannot be taken up. In order to avoid thermodynamically infeasible ATP synthesis the ATP hydrolysis reaction was non-negative by design. The fluxes in Figure 4 for glutamine and CO 2 can be both negative and positive due to these compounds being present in the medium. If the lower end of its bar in Figure 4 is positive that compound must be produced for the cell to grow at maximal growth rate: it is a primary metabolite. When the upper end of the bar is negative it indicates that the compound needs to be taken up for maximal growth.
In Figure 4 we see that in the media where glucose is present (M1-M4) some glucose uptake is essential for attaining the maximal growth rate. In M5-M6 glucose uptake is always zero due to its absence from the medium. Phenylalanine, an essential amino acid, functions as a positive control. Its uptake proves indeed essential in all media for both cell lines as expected. In most media its uptake rate can vary from the amount of phenylalanine in biomass to 0.4, the rate at which it can be used to provide nitrogen to other parts of anabolism. This maximum uptake rate of 0.4 is the same for all media and both cell lines, reflecting that this corresponds to its concentration in the media. In the absence of glutamine and in the presence of low glucose (M4), PLC needs to make full use of this phenylalanine in order to achieve its maximum growth rate, but Huh7 cells could still vary the amount of phenylalanine used whilst attaining the same growth rate. Maximum lactate, pyruvate and ATP production capabilities track the total amount of carbon in the medium within each cell line, with subtle differences between the two cell lines. The gene expression levels appear to be consistent with shifting to virtually complete metabolism of glucose and glutamine to lactate whilst maintaining the maximum growth rate. At the same biomass production flux, lactate efflux could also be as high as 55 mM, corresponding to 2 lactate per maximum glucose consumed plus one lactate per maximal glutamine consumed. However, in all models, the lactate secretion can also be zero while maintaining maximal growth. This shows that glucose conversion to lactate can vary greatly, and may also reflect that in our models the cells can produce and secrete other compounds such as pyruvate. The in silico cells are not addicted to the Warburg effect.
Glutamine uptake is only essential for both PLC and Huh7 in medium 5 and for PLC only in medium 3. In medium 3 for PLC it is then required at a very low amount to achieve optimal growth (as indicated by the upper end of the bar) whereas in medium 5 in both cell lines the maximal uptake bound has to be hit to achieve maximal growth (see the markers at -4). Because (in silico) the maximum growth rate of Huh7 is lower it has the luxury of producing glutamine from glucose whilst growing maximally in M3 whereas this is not possible for PLC. Both cell lines may produce glutamine in M1 and M2 owing to the excess glucose in those media. We conclude that the model cells are insensitive to glutamine concentrations in the medium in the presence of high glucose but glutamine-sensitive in the absence of glucose.
In all models a small amount of oxygen must be taken up. CO 2 (either as CO 2 or as bicarbonate) may either be produced or taken up in M1-M3 for both cell lines and M4 for Huh7 and may only be produced in M4 for PLC and M5-M6 for both cell lines. CO 2 uptake might have to do with the reversal of the isocitrate dehydrogenase reaction, which produces isocitrate as substrate for ATP citrate lyase producing cytosolic acetyl CoA for lipid and cholesterol synthesis [27]. For conditions where CO 2 or bicarbonate production is required this may point to oxidative phosphorylation being required for maximal biomass production.
Oxygen uptake was essential even in conditions where oxidative phosphorylation (interpreted as CO 2 production) seems not to be required. In these cases, oxygen uptake may be necessary for the synthesis of tyrosine, or of cholesterol and other lipids that are part of the biomass definition and absent from our growth media. We checked that removing cholesterol from the biomass equation reduced the need for oxygen, but it did not remove it.
The possibility to grow at maximum rate in the absence of CO 2 production in some conditions, highlights the possibility for cells to obtain all the Gibbs energy they need for maximal growth only from the conversion of glucose to lactate. This may underlie the selection of the Warburg effect by a-social cells [27]. It does not quite correspond to the Warburg and WarburQ effects however: the in silico cells are not addicted to the absence of respiration, as they can still respire all this substrate whilst growing at the same rate. In PLC cells, but not in Huh7 cells, the maximal growth rate in low glucose medium without glutamine (M4) does require oxidative phosphorylation, consistent with the glutamine to lactate pathway elucidated by Damiani et al [27]. These and other apparently minor differences between cell lines in our FVA results are of interest, as they suggest that drugs, in this case ones that inhibit respiration, should be effective against some cancer cells and not others, also depending on extracellular metabolic conditions.  Table 1). In these calculations ATP is treated differently from the others. For the others the reaction was and remained present and potentially carrying varying flux when any of the yet other fluxes was manipulated in the FVA. For the ATP, the ATPase reaction was absent (no growth rate independent maintenance therefore) when doing FVA for any of the others and only present when ATP synthesis was manipulated by forcing flux through an ATP hydrolyzing added reaction. For visibility we used starred markers to indicate flux ranges with a width less than 0.1. These markers are typically located at minimal or maximal flux boundaries: e.g. zero glucose uptake in M5 and M6 or maximal phenylalanine uptake in PLC M4.
We further explored this by plotting some essential reactions for respiration to occur analogously to Figure 4 in terms of their minimal and maximal possible flux allowed while maintaining maximal biomass flux for the medium and cell type specified ( Figure 5). In media with glucose the maximum growth rate does not require flux through cytochrome oxidase and oxidative phosphorylation with the exception of M4 for PLC. Figure 5 suggests that for PLC respiration in terms of flux through cytochrome oxidase is required for maintaining the maximal biomass flux in media 4-6 whereas for Huh7 this is required for media 5 and 6: rather than a range of fluxes, a precise non-zero flux magnitude is required. This suggest that only in those cases of limiting metabolic substrate, the maximal growth rate depends strictly on ATP produced by oxidative phosphorylation. This is in full agreement with the interpretation of the CO 2 + bicarbonate panel in Figure 4. In all other cases respiration is optional for maximum biomass synthesis flux, suggesting that the cells can obtain their ATP from other processes including aerobic glycolysis. To maintain their maximum growth rate at the 5.6 mM glucose concentration, they do need to use virtually all that glucose however ( Figure 4).
Because of its assumptions of maintenance of maximum growth rate and full capability of the network to allow for various fluxes, flux variability analysis makes few predictions that may be put to the test in this study. Exceptions are (i) both cell lines should be capable of consuming the 5.6 or 25 mM glucose offered to them, (ii) they are not addicted to a 100% aerobic glycolysis, but can reduce lactate production without giving up their maximum growth rate, (iii) at glucose concentrations around 5 mM they would make use of all that glucose to grow maximally, (iv) They should be capable of catabolizing glutamine both in the absence and presence of glucose.  Table S3 for the detailed reactions.
Uptake and secretion fluxes of metabolites have units of mM/h. Since we set maximal fluxes equal to medium concentrations (Table 1) these rates differ from reality by some factor.

Growth rate
In a first attempt to examine the applicability of our computational results, and to explore experimentally whether the concentrations of glucose and glutamine in the media control the cells proliferation rate, we evaluated the doubling time for Huh7 and PLC cell lines in six media containing different amounts of glucose and glutamine (see Figure 6A the Huh7 cell line. In both HCC cell lines, the specific growth rate was almost the same for low and high glucose levels, 5.6 and 25 mM respectively. In the absence of glucose, cells grew more than 5 times more slowly. In the presence of glucose, the growth rate did not correlate with the level of glutamine in the medium, suggesting a lack of glutamine addiction in both cell lines. Only for media without glucose, the addition of glutamine to the growth medium led to a higher than 50% increase in the growth rate of both cell lines.

Glucose and glutamine consumption
During the same first 24 hours as when we measured growth rates, we observed a decrease in the concentration of glucose with time such that some 5 mM was consumed by the PLC cells ( Figure 6B). In the PLC cultures that started with 5.6 mM glucose this resulted in almost full glucose depletion after the first day and night. In case of the Huh7 cell line, the consumption of glucose was higher in rich glucose medium than in low glucose medium reaching some 8.5 mM and 5 mM, respectively. Further velocity of glucose consumption for Huh7 was the same, while in case of PLC cells, after 24 h of culturing, stabilization appeared and lasted for next 20 hours ( Figure 6B and 6C). We did not observe a difference in glucose consumption between media with and without glutamine. With respect to glutamine consumption by the cells we found that the level of glutamine did not change maintaining the 4.5 mM level for media with glutamine (that contained originally 4 mM) and 0.5 mM without. i.e. 35% (Figure 7). For the Huh7 cells the trends were similar in case of M1 -M3 medium.
In case of medium M4 such production was achieved only after 45 hours. In M1 and M2 lactate increased linearly with time attaining 10 mM at the end of experiment. Consistently, we did not observe any difference between M1 and M2 for both cell lines nor did we observe such differences between M3 and M4 for PLC; for Huh7 small deviations were noticed. For the two cell lines we examined, this proves an independence of lactate production from glutamine access. percentage of glucose converted into and secreted as lactate is substantial even when there is sufficient oxygen for the glucose to be completely oxidized to carbon dioxide. For the maximal aerobic glucose fermentation, the value of P should be close to 100%, only deviating from that number because of carbon ending up in daughter cells. Values higher than the 100% minus the correction for growth rate, or values higher in glutamine containing media than in media without glutamine could reflect lactate production from glutamine, i.e.
the WarburQ effect.

Discussion
Many human cancers including HCC produce lactate aerobically. This often correlates with tumor progression and worse clinical outcomes in cancer patients. Understanding how the Warburg effect is regulated in cancer and what nutrients are essential for tumor cells to survive, grow and divide may be relevant for identifying new therapeutic interventions. In this work we investigated the effect of two most suspected nutrients i.e. glucose and glutamine by way of six different medium conditions (see Table 1 and 2) using both computational and experimental approaches.
In order to perform the in silico analysis we elaborated a new method to generate models based on the latest human metabolic map (i.e. Recon3D [37]) that take transcriptomic data as well as different medium conditions into account. We thereby obtained twelve (six for PLC and six for Huh7) models that integrate transcriptomic (nature) and nutrient concentration (nurture) data through a new FBA approach, by using proposed by us a new NNS factor. Our new methodology enabled us to predict metabolic changes upon nutrition shift by flux balance analysis.
In this work, we studied the Warburg effect per se. We therefore selected cell types that do exhibit the Warburg effect clearly. We used transcriptomic data from the cell population in its entirely and the models are a representation thereof. Single cell transcriptomic studies suggest that cancer cell populations are heterogeneous, but single cell metabolomics does not yet enable us to examine the consequences for metabolism.
Thus, we assumed in model that the cell population was homogeneous.
Our findings offer support for the predictive potential of genome-scale metabolic maps together with transcriptomic data sets and nutrients composition. The predictions address the carbon and energy metabolism of cancer cells. Because much of this metabolism is essential for cell survival, the potential may translate to new drug targets in a long-neglected area of drug discovery. Indeed, we have shown that dependencies on medium composition and gene expression can be computed and that the results are then relevant enough to be compared with experimental work. Moreover, in most cases the predictions were matched by the experimental results.
Notwithstanding these successes, our methodology comes with a number of issues.
One of these relates to whether expression level information was properly translated to flux.
In contrast to several methods developed to extract context-specific models [60][61][62] that focused on threshold selection with exception of some essential metabolic functions that are needed for cell growth, our method used a linear relationship between flux bounds and the transcriptome. Such a linear approach was first proposed by Colijn et al. [55] and applied in the MaREA framework by Graudenzi et al. [54], where it was shown to be of use in comparing the metabolism of samples in distinct subgroups. The approach assumes that the V max of a reaction is proportional to the level of the mRNA encoding the enzyme, with a proportionality constant equal for all reactions. It thereby neglects differential translation and posttranslational regulation, and assumes all k cat to be equal. Both assumptions are often far from the truth but have the advantage that they can be improved upon in the future by consulting more literature information for the steps that turn out to be important. Keeping the transcriptomics data the same as in the literature source, we were further assuming that the transcriptome remained constant at the time scale of 12 hours that essentially determined our results; we propose that a single cell cycle or less (i.e. < 24 hours) should qualify for this assumption. Furthermore, we assume that the cell lines in our experiments have not significantly evolved compared to those used in the transcriptomics datasets [46].
We did not take into consideration epigenetic and gene-expression changes that could occur  showed that in the absence of glucose, the cells should be able to grow on glutamine but with a much lower growth rate than on glucose alone.
According to the flux variability analysis predictions both cell lines should have a higher potential; they should be capable of consuming the 5.6 or 25 mM glucose offered to them. Yet, our in vitro observations indicated that PLC cells confronted with the 25 mM did not consume more than the around 5 mM during the first 24 hours. They did not make use of that full potential. We conclude that there is a maximum amount of glucose per unit time that the cells 'wish to' handle at glucose levels above a few mM, because other issues than energy and carbon may limit the cells' growth. We use the term 'wish to' to indicate that this may be an issue of metabolic regulation: the gene expression would allow for higher glucose uptake fluxes.
The cells that had consumed in the first 24 h some 5 mM of the 5.6 mM glucose they had been incubated with, continued to grow for the next 24 hours at virtually the same growth rate; they must have done this at a much-reduced glucose consumption rate. Since they also stopped producing lactate, we suspect that the small amount (approximately 0.5 mM) of glucose left provided the cells with sufficient ATP to drive their continued growth. The cells may have reverted from lactate production to respiration with its more than 15-fold higher ATP yield. This highlights that there may be a limitation to these cells' addiction to lactate production: these cancer cells can shift to glucose oxidation. We observed a slightly different effect in the case of the Huh7 cell line: the consumption of glucose was almost two times higher in rich glucose medium (25 mM) than low on (5 mM) during the first 24 hours.
In silico analyses predicted and in vitro experiments confirmed that both cell lines, at glucose concentrations around 5 mM, made use of all available glucose to grow maximally. They also showed that upon glucose depletion and if asked to grow at maximal growth rate, the cells should shift to glucose respiration (Figure 4 and 5).
Lactate production correlated with glucose consumption. Between 15% and 44% of the glucose was transformed to lactate during first 24 hours. Partially contrary to the predictions by our Flux Balance Analysis, we did not detect consumption of glutamine in any medium and for either cell line in the experimental work: the glutamine concentration was constant through all time points, maintaining the 4.5 mM level for media with extra glutamine (4 mM) and 0.5 mM for the media without, even though in absence of glucose from the medium, the cells grew faster when they did have access to glutamine. This observation was in contrast with the modelling results according to which, the cell lines should be capable of catabolizing glutamine both in the absence and presence of glucose. According to the computational experiments, glutamine consumption should have been required to achieve optimal growth in medium 3 for PLC and in medium 5 for both PLC and Huh7 cells. Here the limitation of FBA that it requires an objective function, such as the assumption that the cells have been reprogrammed to provide a maximal growth rate, may be at fault.
Under aerobic conditions, 'social' cells use mitochondria to oxidize the glycolysisderived pyruvate. In the absence of oxygen however, they ferment glucose into lactate. The phenomenon of aerobic lactate production from glucose has been termed the Warburg effect and is characteristic of many fast-growing cancer cells [10,12, [62][63][64][65]. DeBerardinis et al. reported that slow-growing rat hepatoma cells were more oxidative, whereas the more proliferative hepatomas were more fermentative [69]. Zu and Guppy reported that ATP derived through aerobic lactate production in various cancers and cell lines accounts for only 17% of the total ATP [70], still more however than the expected 5% expected for equal lactate production and respiration.
While the involvement of glucose in cancer cell metabolism has been widely studied [10,62,63,71], most such cancer cells were grown on cell culture media containing high levels of glucose, i.e. between 10 and 25 mM [72], and without addressing glutamine utilization. Many cancer cell lines, including some of the Warburg type, display addiction to glutamine [23,27] with some dying rapidly if they are deprived of glutamine. But there are also cancer cells that are not addicted to glutamine [66].
The next question that remains to be addressed relates to the heterogeneity of population and the Intracellular Lactate Shuttle [10b,12b]. Future investigations with single cell transcriptomic/metabolomic data and multi-cells genome scale model are warranted to address this question. In summary, our novel approach that integrates nature (the expression level) and nurture (the environment), let us to corelate proliferation, glucose consumption, and lactate production with the medium concentrations of glucose and glutamine. Our findings delineate a limitation of glucose consumption by HCC cells. Their growth rate was completely independent of glutamine in case of media containing glucose. In view of the requirement of nitrogen for biomass synthesis, this reflects the importance of gratuitous components of growth media. We show that even under typical physiologically relevant conditions, where plasma glucose lies between 4 and 6 mM, production of lactate occurs pointing to aerobic glycolysis. Our results suggest that only a reduction of blood glucose levels to below 5 mM may result in decreased cancer cell proliferation. This conclusion is supported by recent research where Metformin (a drug that decreases blood glucose levels mainly by suppressing liver glucose production through gluconeogenesis) turned out to be more effective in enhancing chemotherapy sensitivity of cancer cells under reduced glucose conditions [72].

Author Contributions
EWT formulated the problem, proposed and designed the project, with advice from HVW.
EWT conceived the strategy of the study with contributions from TDM and HVW. EWT designed in vitro experimental work that was performed by EWT and DP. Computational analyses were designed and performed by TDM with advice from EWT and HVW and interpreted by TDM, EWT and HVW. EWT, TDM and HVW wrote the manuscript.

Conflicts of Interest
The authors declare that there is no conflict of interest.
with early modeling work. We thank Chiara Damiani for discussing possible extensions of the MaREA methodology, inspiring us to the implementation here.

Supplementary Files
• Supplementary material • Supplementary Excel