Multi-scale, whole-system models of liver metabolic adaptation to fat and sugar in non-alcoholic fatty liver disease

Non-alcoholic fatty liver disease (NAFLD) is a serious public health issue associated with high fat, high sugar diets. However, the molecular mechanisms mediating NAFLD pathogenesis are only partially understood. Here we adopt an iterative multi-scale, systems biology approach coupled to in vitro experimentation to investigate the roles of sugar and fat metabolism in NAFLD pathogenesis. The use of fructose as a sweetening agent is controversial; to explore this, we developed a predictive model of human monosaccharide transport, signalling and metabolism. The resulting quantitative model comprising a kinetic model describing monosaccharide transport and insulin signalling integrated with a hepatocyte-specific genome-scale metabolic network (GSMN). Differential kinetics for the utilisation of glucose and fructose were predicted, but the resultant triacylglycerol production was predicted to be similar for monosaccharides; these predictions were verified by in vitro data. The role of physiological adaptation to lipid overload was explored through the comprehensive reconstruction of the peroxisome proliferator activated receptor alpha (PPARα) regulome integrated with a hepatocyte-specific GSMN. The resulting qualitative model reproduced metabolic responses to increased fatty acid levels and mimicked lipid loading in vitro. The model predicted that activation of PPARα by lipids produces a biphasic response, which initially exacerbates steatosis. Our data support the evidence that it is the quantity of sugar rather than the type that is critical in driving the steatotic response. Furthermore, we predict PPARα-mediated adaptations to hepatic lipid overload, shedding light on potential challenges for the use of PPARα agonists to treat NAFLD.


SUPPLEMENTARY INFORMATION
Multi-scale, whole-system models of liver metabolism adaptation to fat and sugar in nonalcoholic fatty liver disease. Maldonado

Multi-Scale Model of Hepatic Monosaccharide Metabolism Reconstruction and QSSPN Simulation
The multi-scale model of the regulation of sugar uptake and its effect on global hepatocyte metabolism is composed of a fully parameterised ordinary differential equation (ODE) model of insulin signalling 1 and a hepatic genome scale metabolic network (GSMN) 2 constrained by consumption and release flux bounds derived from a metabolomics in vitro study. 3 A Petri net formalism was used to represent the signalling network composed of monosaccharide transport and insulin signalling (Fig. S2). Reconstruction was performed with the Petri net editor software Snoopy 2 and run on Mac OS X; v1. 13 were also used to link qualitative and quantitative aspects of the DT (e.g. G6Pase). Transitions can be used to define the reaction rate constant for a reaction through the 'RATE' comment, with a rate of 1 being used as default if this is not present. In addition, the 'FLUX' comment can be used to extract real flux values from the simulated solution; 6 these transitions are indicated as red rectangles within Figure S2. Standard (point-ended arrows) and read (circleended arrows) edges were used to simulate consumption or catalyzation of the interactions, respectively. Edges were tagged to specify the activity of the pre-place, e.g. 'ACTIVITY 0' followed by 'END' denotes mass action activity. To represent in vivo-like transportation of glucose and fructose, Michaelis-Menten kinetics (Eq. 1) were utilised to calculate the fluxes within the activity list of the monosaccharide constraint places, such that Eq. 1 where is the rate of the reaction, is the maximum rate possible for the reaction, [ ] is the substrate concentration and KM is the Michaelis-Menten constant which is unique to each substrate.
The was set at 5 mmol/g DW/h for both glucose and fructose based on favouring whole cell uptake, rather than recombinant systems. 7,8,9,10 HepG2 cells are known to express more than one glucose transporter, e.g. GLUT1, GLUT2, GLUT3 and GLUT9, 11,12,13,14 but not the main fructose transporter, GLUT5. 12,15 Glucose transport is symmetrical. 9 Therefore, the activity list for glucose and fructose transport constraint places are symmetrical in regards to the maximal import and export rates represented as lower and upper bounds, respectively.
Monosaccharide transport rates 7,8,10 were converted to mmol/g DW/h based on assumptions outlined in Table S1 and  . S3B) resulted in similar behaviour as seen in the original publication (Fig. S3A). The insulin model was converted into a Petri net formalism using Snoopy (Fig. S2F). Large places denoted signalling molecules, medium sized places represent reaction rates and the smallest places represent degradation. Places that are shaded grey represent clones of individual places to enhance visualisation: all clones for a single entity contain the same concentration/flux. Transitions were labelled as the reaction rate IDs in the original publication, and were parametrised by their pre-place. The insulin model simulated in QSSPN was consistent with the model simulated in COPASI, producing identical results ( Fig. S3B, C).
The kinetic insulin model was then coupled with the HepatoNet1 GSMN. Within the original kinetic insulin model, only one gene product had a direct link to a metabolic reaction in HepatoNet1: G6Pase (r0396, H2O(r) + Glucose-6P(r) → Glucose(r) + Pi(r)). To enhance the regulatory connection between the insulin regulatory network and the metabolic network, phosphorylated glycogen synthase kinase (pGSK)3β was designated as the active glycogen synthase (GYS; r1388, 3 UDP-glucose(c) + Glycogenin-G8(c) → Glycogenin-G11(c) + 3 UDP(c)). 20,21 Comments in the original publication also noted that another compound, phosphoenolpyruvate carboxykinase (PEPCK), had a similar behaviour to G6Pase but it was not included in the original kinetic model. 21,22 Therefore, the activity of G6Pase was used to represent the activity of PEPCK and connected with HepatoNet1 reactions r0123-4 for both cytosolic (c) and mitochondrial (m) compartments of reaction; GTP + OAA → GDP + PEP + CO2.
The HepatoNet1 GSMN, 2 representing human liver whole-cell metabolism, was downloaded from the BioModels database (MODEL1009150000, SBML L2 V4). 23 All fluxes were checked to ensure the direction of reactions in the stoichiometric matrix. All reversible fluxes were adjusted to a minimum flux of -1 and a maximum flux of 1; irreversible reaction fluxes were adjusted to 0 to 1 or -1 to 0 depending on directionality.
Specific adjustments were made to HepatoNet1 to better detail metabolic pathways of interest accurately. The polyol pathway was included as this has been found to be a potential influence on the development of hepatic steatosis 24 and is part of fructose metabolism in the liver. 25  also added to represent the hepatic fructose sodium co-transport reaction. 27 Thus, the total number of metabolites is 778 and the total number of reactions is 2542 within the modified GSMN.
A biomass function was used as a simulation constraint to ensure the presence of the essential building blocks (amino acids, nucleoside triphosphates and deoxynucleotide triphosphates) for biomass were available; representing the basic metabolic requirements of a human cell. This function was adapted based on a published whole human metabolism GSMN, Recon 2, biomass 28 (Table S5). FBA simulations before and after the modifications listed above, and using the same external metabolite constraints (Table S3), were performed: BIOMASS optimal maximisation flux was not altered by the modifications: FBA = 0.03152.
Biomass doubling time estimations for HepG2 cells range from approximately 20 to 60 hours. 29,30,31,32 To reflect this, the BIOMASS lower and upper flux bounds were set at 0.01666667 and 1000.0, respectively, representing a minimal doubling time of 60 hours.
Additional constraints were included to represent the physiologically relevant kinetic activity of the first steps of glucose and fructose metabolism. Flux constraints (mmol/g DW/hour) were set as a ratio of fluxes. As a conservative measure, flux ratio assumptions were based on a hexokinase Michaelis-Menten constant KM for glucose and fructose phosphorylation at carbon-6, as hexokinase II is the most likely expressed hexokinase in HepG2 cells. 14 The hexokinase KM for glucose and fructose phosphorylation were estimated as 0.047 mM and 11.5 mM, respectively. 33 Thus, hexokinase II flux ratio is set for the phosphorylation of glucose specific reactions (r0353-5) at 1, and reactions specific to fructose phosphorylation (r0356-8) was set at 0.005. The polyol pathway flux was also adjusted. The KM for glucose phosphorylation by hexokinase was estimated to be 0.047 mM, while the KM of the first reaction in the polyol pathway by aldose reductase enzyme was estimated as 65.8 mM. 34,35,36 Thus, hexokinase phosphorylation of glucose was set as 1 (as above) and the flux of the reaction specific to AKR1 was set as 7.14 x 10 -4 . For computational purposes, the TAG production reaction (r1223) was adjusted as a positive flux as QSSPN simulations can only maximise positive objective function flux values. Additionally, TAG production reaction (r1223) lower and upper bound values were set as 0 and 10, respectively, so as not to mask the optimisation flux if above 1.
To represent the nutrient composition of the culture medium, transport reactions were included as an external exchange set. The HepatoNet1 physiological import and physiological export set (PIPES) was modified to represent the nutrient composition of cell culture medium as shown in Table S4. The fluxes of the external metabolite set were constrained by setting the lower and upper bound fluxes as the maximum consumption and release rates derived from the NCI-60 cell lines consumption and release of metabolites dataset (CORE) 3 , respectively, as previously done in 37 .

Model Simulation
The QSSPN algorithm integrates QSSF set by the metabolic network with DT represented by the Petri net network. 5 vi.
The majority of dFVA simulations were run with a maximal time step of 0.01. However, this proved to be computationally expensive when simulating the integrated model of insulin and HepatoNet1 optimising for TAG production (as the objective function). Therefore, this setting was adjusted to 0.1. Under these conditions, simulations were still of sufficiently high resolution to allow for examination of the network behaviour.

Sugar Consumption Assay
Sugar consumption was monitored by sampling culture media over time and quantified by using glucose and fructose assay kits (Abcam, UK). Both glucose and fructose were stable in medium over the experimental period, and no extracellular degradation was expected. The disappearance of monosaccharide from the medium is highly suggestive and was regarded as cellular uptake and consumption. Briefly, medium or diluted medium (medium:assay buffer,

Pathway Enrichment Analysis
Enriched pathways (P<0.05), within the transcriptomics and proteomics datasets described below, were identified against the KEGG 39 and BIOCARTA 40 databases using the DAVID 41 online bioinformatics suite. Results from these analyses were systematically combined and quantitatively represented using a hive plot formalism. 42

Transcriptomics
Microarray experiments were performed as previously described. 43 Briefly, Huh7 cells were incubated with 100 µM or 200 µM of palmitic acid in serum-free DMEM for 24 hours. Cells were homogenised in Trizol (ThermoFisher, UK) and total RNA isolated as per the manufacturer's instructions. Total RNA quantity and purity was assessed by absorbance spectroscopy (λ 260nm, 280nm and 230nm); integrity was confirmed using the Agilent 2100 bioanalyser and the RNA 600 LabChip kit prior to hybridisation to human gene expression microarrays (Agilent Technologies, UK). Differential expression of genes (≥1.5 fold relative to controls, P<0.05) were identified in both treatment groups and taken forward for pathway enrichment analysis.

Proteomics
Details of the in vivo proteomics study have been described in detail elsewhere. 44 Briefly, C57BL6 and APOE-/-mice were fed either high-fat or normal chow diets for 12 weeks after which animals were sacrificed, livers were resected, homogenised and protein isolated and enriched for cytoplasmic and membrane proteins. Fractionated samples were iTRAQ labelled for proteomics analysis by mass spectrometry. Proteins identified as being differentially expressed (relative to wild-type animals fed a normal diet) were carried forward for pathway enrichment analysis.

PPARα Model Simulation
Simulation of the PPARa model was performed as previously published 5 with the exception for updates to the gene expression model (Fig. S4), the biomass function (Table S5), and exchange set (Table S4) detailed herein.

SUPPLEMENTARY FIGURES FIGURE S1
Figure S1. Modelling strategy overview. The computational model construction incorporates information from both in-house experimental data and literature. Although not exclusively, literature-based data was incorporated into our quantitative kinetic regulatory network, while in-house omics data was also used to aid the qualitative gene regulatory network construction.
The quantitative and qualitative model is simulated to deliver outputs that can be used to compare to and/or inform experimental results. This is an iterative modelling approach as experimental results can then be used to update model design to deliver predictions more representative of known biology.  Table S3. B A biomass constraint to represent a 'healthy hepatocyte' (Table  S4). C Monitoring outputs, e.g. TAG production. D The objective function. E The sugar transport system. F The insulin kinetic model reconstructed in a Petri net formalism with coloured ovals used to highlight modules used by Kubota and colleagues (2012) 1 ; the PEPCK module in black oval was interpreted as described in supplementary methods.  HepG2 viability by LDH (n=2) and MTT (n=8); HuH7 by LDH in low (n=3) and high (n=4) glucose; HuH7 by MTT in low (n=4) and high (n=5) glucose. Data are relative to serum cultured cells set at 100 % (dotted line) and shown as mean ± SEM. Two-way ANOVA with Bonferroni's test post hoc was used to detect differences between serum-free and added 2% DMSO treated cells with statistical significance indicated by * P < 0.05 and **** P < 0.0001.      3 . A negative flux represents import while a positive flux represents export, and are set as the lower and upper bounds, respectively. Metabolite fluxes that were not measured by 3 were set accordingly as described by HepatoNet1 PIPES. 1 Metabolites were not measured in 3