Formal modeling and analysis of the hexosamine biosynthetic pathway: role of O-linked N-acetylglucosamine transferase in oncogenesis and cancer progression

Metabolic reprogramming and oncogenesis

The capability of tumors to reprogram cellular metabolism and promote uncontrolled proliferation has only recently emerged as a cancer hallmark (Hanahan & Weinberg, 2011; Cairns, Harris & Mak, 2011; Dang, 2012). However, alterations of energy metabolism in cancer cells to stimulate cell growth and division first observed by Otto Warburg date back to the early twentieth century (Warburg, 1910; Weinhouse et al., 1956; Warburg, 1956). To fuel growth, tumor cells flip a metabolic switch to reprogram glucose metabolism from oxidative phosphorylation to aerobic glycolysis, and secrete lactate (“Warburg-effect”). While glycolysis produces adenosine 5′-triphosphate (ATP) faster, this process occurs at a less efficient rate. To compensate for this metabolic switch, neoplastic cells become addicted to glucose and glutamine to maintain rapid cell proliferation (Potter, 1958; Vander Heiden, Cantley & Thompson, 2009; DeBerardinis et al., 2007; Lunt & Vander Heiden, 2011; Mullen et al., 2012; DeBerardinis et al., 2008). This phenomenon of increased glycolytic flux and glucose uptake importantly increases flux into biosynthetic pathways, such as the Hexosamine Biosynthetic Pathway (HBP) (Fig. 1).

Studies have also found an association between glycolytic fueling and mutant tumor suppressors or activated oncogenes that play critical roles in evading apoptosis and promoting proliferation of tumor cells (Jang, Kim & Lee, 2013; Ying et al., 2012; Gross, Van den Heuvel & Birnbaum, 2008; Haq et al., 2013). Ying et al. (2012) found that over-expression of Kras was implicated in tumor initiation by controlling tumor metabolism and channeling glucose intermediates into the HBP. Kras upregulation of glycolytic enzymes, glucose transporters, and glutamine: fructose-6-phosphate amidotransferase 1 (GFAT1), drives increased HBP flux and cellular O-linked N-acetylglucosaminylation glycosylation (O-GlcNAcylation) in cancer (Hsu & Sabatini, 2008; Ma & Vosseller, 2013). A fundamental role of the HBP is to control O-GlcNAcylation. O-GlcNAcylation is a post-translational modification catalyzed by O-GlcNAc-transferase (OGT) that covalently adds a GlcNAc sugar moiety to hydroxyl groups of serine/threonine residues of proteins (Torres & Hart, 1984; Hart, Housley & Slawson, 2007; Issad & Kuo, 2008). This process is antagonized by O-GlcNAcase (OGA), which allows for dynamic regulation of O-GlcNAcylation in cells (Fig. 1).

Recently, O-GlcNAcylation has been proposed as a novel cancer hallmark and approach for cancer treatment due to its significant regulatory role in tumorigenesis (Fardini et al., 2013). Increased O-GlcNAcylation, termed hyper O-GlcNAcylation, and elevated OGT levels have been observed in various tumor types, including cancers of the breast, lung, liver, bladder, endometrial, prostate, pancreas, and colon (Ying et al., 2012; Gu et al., 2010; Mi et al., 2011; Zhu et al., 2012; Rozanski et al., 2012; Krześlak et al., 2012b; Lynch et al., 2012). Importantly, the inhibition of OGT has been associated with decreased proliferation of breast and prostate cancer cells (Caldwell et al., 2010; Itkonen et al., 2013). Understanding the proteins that control deregulation of cellular energy metabolism and hyper O-GlcNAcylation is needed to elucidate the mechanisms of metabolic switch in cancer cells, characterize the glycolytic phenotype, and decipher the link to cellular growth and apoptotic pathways.

Signaling pathways

The ability of tumors to promote an inflammatory response and escape immune destruction also enables cellular proliferation and evasion of innate immunity (Dvorak, 1986; Colotta et al., 2009; Hanahan & Weinberg, 2011; Markert, Levine & Vazquez, 2012; Kroemer & Pouyssegur, 2008). Inflammation is the protective response of the innate immune system to a physiological, physical, and/or oxidative stress. The development of innate immunity is associated with the NF-κB signaling cascade, where NF-κB is activated through subunits of the IKK complex in response to stimuli (Karin, 2009; Hoesel & Schmid, 2013). Activation of NF-κB targets and eliminates transformed cells, Disis (2010) supporting subsequent increases in apoptotic processes as an inflammatory response (Ernst, 1999; Cordon-Cardo & Prives, 1999). Further studies revealed additional roles of NF-κB in controlling normal cellular and malignant processes, such as proliferation, apoptosis, and metabolism (Guttridge et al., 1999; La Rosa, Pierce & Sonenshein, 1994; Perkins, 1997; Moretti et al., 2012; Kawauchi et al., 2009; Kawauchi et al., 2008).

The subunits of NF-κB contain sites for post-translational modifications to promote cross-talk with signaling pathways. O-GlcNAcylation of the c-Rel subunit of NF-κB was recently demonstrated to mediate the expression of various cytokine-encoding genes involved in the immune response (Alexandrov et al., 2013). More recent studies have noted that expression of OGT, the enzyme that catalyzes O-GlcNAcylation, is correlated with c-Myc protein levels and may be involved in protein stabilization (Itkonen et al., 2013).

c-Myc belongs to the PI3K-mTOR-MYC signaling pathway, one of the most commonly mutated pathways in cancer (Fig. 1). Activation of this signaling cascade has been shown to increase hyper O-GlcNAcylation activity in breast cancer (Sodi et al., 2015). Notably, treatment of tumor cells with PI3K and mTOR inhibitors led to decreased protein expression of OGT and overall lower levels of O-GlcNAcylation. The PI3K pathway has also been reported to cross-talk with Forkhead box M1 (FoxM1), an oncogenic transcription factor (Major, Lepe & Costa, 2004). FoxM1 plays a critical role in cancer metabolism, as the reduction of O-GlcNAc levels and OGT in cancer cells is associated with a decrease in protein expression of FoxM1 in breast cancer (Caldwell et al., 2010). However, studies to date have not detected O-GlcNAc modifications on FoxM1.

In addition to NF-κB and c-Myc, p53 has also been shown to be directly O-GlcNAcylated. p53 plays critical roles in DNA damage repair and apoptosis, and is one of the most frequently mutated genes in cancer. Hyper O-GlcNAcylation of p53 stabilizes the tumor suppressor and decreases p53-MDM2 interaction to block proteolysis (Fig. 1) (Yang et al., 2006). In contrast, overexpression of OGA, the antagonist to OGT, stimulates MDM2-p300 interaction and degrades p53 (Soesanto et al., 2008). Consequently, loss of p53 activates NF-κB to increase aerobic glycolysis and support tumor metabolism (Kawauchi et al., 2008). Under stress, p53 can induce cyclin-dependent kinase inhibitor p21 to arrest the cell cycle and inhibit proliferation (Gartel & Tyner, 1999). Taken together, disentangling the complex interplay between NF-κB, c-Myc, p53, MDM2, FoxM1, p21, and OGT is critical to understanding the roles of hyper O-GlcNAcylation, pathway signaling and cross-talk, metabolism, and programmed cell death in cancer.

Our contribution

In this study, we examine the role of hyper O-GlcNAcylation in cancer progression by regulating the activation of oncogenes. We construct a qualitative Biological Regulatory Network (BRN) comprised of important entities involved in O-GlcNAc signaling to demonstrate activation and inhibition relationships. Model parameters are computed from known experimental observations by using a formal verification technique, called model checking. These parameters are used to translate BRN into a qualitative model which highlights important behaviors as trajectories, stable states, and cycles. Network analysis of the qualitative model is performed to identify important trajectories involved in oncogenic activation, cancer progression, and recovery. We identify significant cycles that represent normal behavior of the overall system and use hybrid modeling to compute delay constraints, which limit the system to maintain homeostasis. A similar modeling approach has been used in the past for formal modeling of biological networks, including the MAL-Associated Biological Regulatory Network (BRN) (Ahmad et al., 2012), the regulatory network of dengue virus pathogenesis and clearance (Aslam et al., 2014), the mechanism of tail resorption in tadpoles (Khalis et al., 2009), and the immunity control mechanism in bacteriophage lambda (Richard, Comet & Bernot, 2006). The results of our study highlight that O-GlcNAc transferase (OGT) plays an important regulatory role in oncogenic activation. The qualitative model reveals that persistent over-expression of OGT and c-Myc leads to deadlock state, from which the system cannot proceed to a recovery state. Another important insight obtained from the model is that silencing of the OGT and c-Myc loop decreases glycolytic flux and results in programmed cell death. Based on the results of network analysis carried out using Cytoscape, we identify a significant cycle, which highlights the important role of p53-Mdm2 oscillations to bring the system towards recovery state. The results of hybrid modeling suggest delay constraints to maintain homeostasis. We compare the important insight gained through computational modeling to show that results are in agreement with previous studies. Together, our findings suggest that the OGT and c-Myc loop is critical in tumor progression, and targeting these mediators may represent a novel therapeutic strategy to regulate hyper-O-GlcNAcylation for the treatment of cancer.

Qualitative modeling framework

The qualitative modeling framework introduced by René Thomas uses a graph-theoratic approach to model BRNs. Each BRN is modeled as a weighted directed graph in which nodes represent biological entities such as genes or proteins, whereas the activation and inhibition relationships between nodes are represented by edges. Here, we briefly introduce semantics of the qualitative modeling framework, mainly adopted from Bernot et al. (Bernot et al., 2007; Bernot et al., 2004; Saeed & Ahmad, 2014).

“

• V is the set of all vertices or nodes

• E is an ordered pair of nodes i.e., if e ∈ E, then e = (v_i, v_j) and $V=\left\{v_i,v_j\right\}$.

The edge (v_i, v_j) is directed from v_i to v_j, where v_i is called the head and v_j is called the tail. In G, the set of predecessors and successors of a node v_j are denoted as $G_{v_j}^{-}$ and $G_{v_j}^{+}$, respectively.

Each node v_i ∈ V has its abstract expression level in the set δ_vi = {0…ℓ_vi} where ${\ell }_{v_i}\le \left|G_{v_i}^{+}\right|$. The state of a BRN is a configuration of expression levels of all entities at a particular time instant.

In a given state, each v_i is regulated by its predecessors $G_{v_i}^{-}$, formally denoted as set of resources, ω_vi, defined as follows:

• ∃ a unique v_xϵV such that $s_{v_x}\ne s_{v_x}^{\prime }$ and $s_{v_x}^{\prime }=s_{v_x}\Rsh K_x\left({\omega }_{v_x}\right)$, and

• ∀ $v_yϵV\setminus \left\{x\right\}{s}_{v_y}^{\prime }=s_{v_y}$.”

Parameter inference using model checking

The dynamics of Thomas’s method are generated by translating the interaction graph to a state transition graph using a set of logical parameters, which are not known in advance. The estimation of model parameters constitutes an important step in qualitative modeling of biological networks. Bernot et al. (2004) introduced a method to decipher these parameters by employing a formal verification approach, called model checking. In this approach, known experimental observations are encoded in a temporal logic framework, called Computation Tree Logic (CTL), and then using the model checker, different parameter combinations are evaluated to finally select parameters which satisfy CTL observations. In CTL, experimental observations are encoded into formulas by using a set of quantifiers which define criteria to explore different states or paths originating from a given state. Here, we provide a brief description of these quantifiers, the detailed semantics of which can be found in Clarke, Grumberg & Peled (1999).

• A: This is a path quantifier which enforces that a given property should hold in all paths originating from the given state. The quantifier itself is read as “For all paths.”

• E: Known as the “Existential Quantifier,” this is also a path quantifier which enforces that a given property must hold in at least one path originating from the given state. The quantifier is read as: “There exists a path.”

• G: This quantifier is known as the “Global Quantifier” and is a state quantifier which enforces that a property holds in all states of a path originating from the given state, inclusive of the given state as well. It is read as: “Globally.”

• F: The “Future Quantifier” is the second state quantifier and enforces that a given property must hold in one of the future states in the path originating from the given state. It is read as: “In future” or “Eventually.” The Future Quantifier also covers the current/given state as well when checking the property.

• X: The “Next Quantifier” is the third state quantifier and enforces that a given property must hold in the immediate successor state. It is read as: “Next.”

SMBioNet (Khalis et al., 2009; Bernot et al., 2004) is a tool for the parameter estimation of biological networks, based on the qualitative formalism of René Thomas (Thomas, 1978; Atkinson, 1965). Given a model of a BRN in the form of Thomas’s network and behavioral properties (observations), expressed as CTL formulas, SMBioNet exhaustively enumerates all compatible parameterizations by generating a state graph for each parameter combination and by verifying the formulas on each state graph. The verification of the CTL property is performed by invoking model checker NuSMV (Cimatti et al., 2002). The parameter combinations are reduced by applying Snoussi and observability constraints (Snoussi & Thomas, 1993). Finally, all the models that satisfy the CTL properties are shortlisted. SMBioNet has been applied in studies such as: tail resorption in tadpole metamorphosis (Khalis et al., 2009), and immunity control in bacteriophage lambda (Richard, Comet & Bernot, 2006).

Network analysis

Graph Theory (Bondy & Murty, 1976) plays an important role in the modeling and analysis of processes in several application areas, including systems biology (Pavlopoulos et al., 2011; Barabasi & Oltvai, 2004; Mason & Verwoerd, 2007). The graph-theoratic approaches are employed to analyze topological and structural parameters of biological networks to discover key properties that provide meaningful insights into the functionality of biological systems. Identification of important nodes in a large biological regulatory network is critical in the understanding of cellular mechanisms. The most widely used measure to compute the ranking of nodes in graph-theoratic models, based on the concept of Centrality (Aittokallio & Schwikowski, 2006; Mason & Verwoerd, 2007), mainly originate from Social Network Analysis (Wasserman, 1994). Centrality Analysis has also been employed to investigate important properties of complex biological regulatory networks (Koschützki & Schreiber, 2008).

Hybrid modeling with delays

Discrete modeling provides useful insights into qualitative dynamics of biological networks. However, an increase or decrease in protein expression, described by a step function, is not coherent with actual changes in protein expression taking place within a cell. The concentration level of a protein, for instance, does not jump from one discrete value to another discrete value. In order to capture the sigmoidal change of protein expressions, Ahmad et al. (2007) introduced a new framework based on piece-wise linear equations. In this framework, states of a system are modeled as discrete locations. Additionally, specialized variables, called clocks, are used to specify constraints for transition from one discrete location to another (Fig. 3). Here, we provide a brief description of the hybrid modeling framework, adopted from Ahmad et al. (2007) and Aslam et al. (2014).

“Clock variables are used to measure the ‘delays’ (the time duration) that needs to pass between two consecutive expression levels. Thus, a clock variable h is associated with each protein in the BRN. The initial values of each h are set to zero, which then approach either d⁺ or d⁻. d⁺ signifies a production delay, that is, the delay required to increase the concentration level of the associated protein by 1. Similarly, d⁻ signifies the degradation delay, that is, the delay to decrease the protein concentration by a single level. The rate of evolution of each h is given by the first order derivative dh∕dt = r where r ∈ {0, 1, − 1} (Ahmad et al., 2007).

In most cases, the exact values of the delays associated with the proteins are not known, which is why unvalued parametric delays are used. Thus, the hybrid model was constructed using the Parametric Bio Linear Hybrid Automaton (Ahmad, 2009) defined below.

Let C⁼(X, P), C^≤(X, P), and C^≥(X, P) be the set of constraints using only =, ≤, and ≥, respectively. Here, X and P are the sets of real valued variables and parameters, respectively.

The Transition System related semantics of the parametric Bio-LHA are given below according to the time domain 𝕋, where 𝕋^∗ = 𝕋∖{0}.

Using the semantics of the Bio-LHA, Ahmad et al. (2007) then defined the temporal state space and the invariance kernel set which have been adapted below.

In the hybrid model of the BRN, we denote ϕ(t) for t ∈ ℝ ≥ 0, while the sequence of points of a trajectory and the set of all points in the state space is denoted by S. A particular trajectory is said to be viable if it remains within a prescribed region known as its viability domain. The state pace is denoted by S. A particular trajectory is said to be viable if it remains within a prescribed region known as its viability domain.

Model construction

From study of the existing literature, we construct a qualitative BRN comprised of nine genes and seventeen interactions (Fig. 4). The BRN is composed of a set of well-known regulatory motifs, each of which give rise to a specific functionality of the system. First, we observe an inhibitory set of genes (OGT vs OGA) that produces a positive feedback loop, known to generate multiple stable states. In practice, a positive feedback circuit is comprised of an even number of negative elements (Plahte, Mestl & Omholt, 1995), and acts as a toggle-switch in which only one of the two genes is expressed at a time (Gardner, Cantor & Collins, 2000). On the other hand, a negative circuit is comprised of an odd number of negative interactions, such as the interaction between p53 and Mdm2, that leads to a periodic behavior or homeostasis. An important oscillatory behavior is characterized by two nested regulatory modules involving PI3K: (i) a positive feedback loop via p21 and NF-κB; and (ii) a positive feedback loop between FoxM1, p21 and NF-κB. The logical analysis of these regulatory motifs provides useful information about the potential behavior of a system. However, functional dynamics of a complex system, that involve both positive and negative circuits, can only be rendered with proper parameter values.

Logical parameters

Table 1 enlists the final values of logical parameters which are used to generate model trajectories as a state transition graph. These parameters are computed from known qualitative observations expressed as CTL formulas (Eqs. (3)–(6)), the most important of which is the change in expression levels of OGT and OGA in different types of cancers, including tumors of the breast and colon (Singh et al., 2015; Fardini et al., 2013; De Queiroz, Carvalho & Dias, 2014). The sub-formulas ψ₁ and ψ₂ represent changes in the HBP pathway, triggered by a change in the expression level of two genes: OGT and OGA. The sub-formula ψ₁ describes behavior of the biological system under enhanced OGT expression leading to a future state in which the expression of oncogenes remain high. The sub-formula ψ₂ describes that there is at-least one trajectory in which expression of oncogenes remain low when OGT is initially not expressed. Finally, the sub-formula ψ₃ represents oscillatory behavior exhibited by the HBP pathway and tumor suppressor proteins p53 and p21. We used SMBioNet software (Bernot et al., 2004; Khalis et al., 2009) to select only those parameters that satisfy the CTL formulas. SMBioNet selected four models (Supplemental Information 2 and 3) which show a single deadlock state (1,0,1,1,1,1,1,0,1) and plausible biological trajectories in cancer progression and recovery. (3)${\displaystyle \begin{array}{c}{\psi }_1=\left(\right.\left(\mathrm{OGT}=1\wedge \mathrm{OGA}=0\right)\Rightarrow \mathbf{EF}\left(\right.\mathbf{AG}\left(\right.\mathrm{OGT}=1\wedge \hfill \\ \mathrm{OGA}=0\wedge \mathrm{PI3K}=1\wedge \mathrm{FOXM}=1\wedge \mathrm{P21}=0\wedge \mathrm{CMyc}=1\left.\right)\left.\right)\left.\right)\hfill \end{array}}$ (4)${\displaystyle \begin{array}{c}{\psi }_2=\left(\right.\left(\mathrm{OGT}=0\wedge \mathrm{OGA}=1\right)\Rightarrow \mathbf{EF}\left(\right.\mathbf{AG}\left(\right.\mathrm{OGT}=0\wedge \hfill \\ \mathrm{OGA}=1\wedge \mathrm{PI3K}=0\wedge \mathrm{FOXM}=0\wedge \mathrm{P21}=1\wedge \mathrm{CMyc}=0\left.\right)\left.\right)\left.\right)\hfill \end{array}}$ (5)${\displaystyle \begin{array}{c}{\psi }_3=\left(\right.\left(\mathrm{OGT}=0\wedge \mathrm{OGA}=1\wedge \mathrm{p53}=1\wedge \mathrm{P21}=1\right)\Rightarrow \hfill \\ \mathbf{EX}\left(\mathbf{EF}\left(\mathrm{OGT}=0\wedge \mathrm{OGA}=1\wedge \mathrm{p53}=1\wedge \mathrm{P21}=1\right)\right)\left.\right)\hfill \end{array}}$ (6)${\displaystyle \begin{array}{c}\mathrm{\Phi }=\left({\psi }_1\right)\wedge \left({\psi }_2\right)\wedge \left({\psi }_3\right).\hfill \end{array}}$

The tendency of each gene to change its expression level is a function of presence or absence of its resources. The change in expression level of a gene can be determined by comparing its current state, at any particular time, with values of logical parameters listed in Table 1. The inferred parameters indicate that in the presence of the PI3K activation signal, NFκB maintains a higher expression level (if already expressed i.e., 1) or shows an increase in its expression level. On the other hand, PI3K shows a rise in expression even in the presence of a p53 inhibition signal when both NFκB and OGT are activating PI3K. The expression level of tumor suppressor protein p53 shows an increase only in the absence of the MDM2 inhibition signal. The inferred parameters show increase in expression level of OGT and C-Myc when there is an activation signal between them. The collective behavior of genes involved in a biological system can only be determined only by analyzing trajectories in a state transition graph.

State transition graph

A state transition graph of the biological regulatory network with 512 nodes and 2,304 edges (Supplemental Information 4) is rendered using Cytoscape software (Shannon et al., 2003) (Fig. 5). The graph is generated from selected logical parameters (Table 1) using GENOTECH software (Ahmad, 2009; Aslam et al., 2014; Ahmad et al., 2012) and states are sorted on the basis of betweenness centrality. A parameterized BRN is also attached for the GINsim Tool software (Supplemental Information 8) (Chaouiya, Naldi & Thieffry, 2012). The deadlock state (1,0,1,1,1,1,1,0,1) (Fig. 5D) shows high expression levels of OGT and oncogenes. The immediate predecessors (up to two levels) of the deadlock state have low betweenness centrality, indicated with circles, having comparatively larger diameters and darker colors. The model also shows several cycles, and rendered as an outermost circle (Fig. 5). These states have high betweenness centrality, represented using circles with smaller diameters and lighter colors. The cycles demonstrate normal behaviors of the system characterized by low expression levels of oncogenes, oscillation of the p53-Mdm2 circuit, and moderate expression levels of tumor suppressor proteins. In the state graph, the state of the system at a particular time is represented by a vector containing expression levels of all entities. The normal state is characterized by low expression levels of OGT, PI3K, and FoxM1 along with the presence of tumor suppressor proteins p53 and p21. This state is represented as a vector (OGT = 0, OGA = 1, PI3K = 0, FoxM1 = 0, p53 = 1, p21 = 1). Conversely, the pathogenic state is characterized by high expression levels of OGT along with PI3K and FoxM1. (OGT = 0, OGA = 1, PI3K = 0, FoxM1 = 0, p53 = 1, p21 = 1) represents a normal or recovery state of the system. Biological systems, under normal circumstances, exhibit oscillatory behavior or homeostasis during which the overall state of the system remains in a cycle of normal states. Therefore, the desirable qualitative model should exhibit pathogenic trajectories along with normal homeostatic behavior represented as a cycle or closed path. The normal behavior, characterized by low expression levels of OGT and the presence of OGA, is encoded as CTL formula ψ₂. It states that, under normal circumstances, when the expression level of OGT is low, the system will always remain in a stable state characterized by low expression of oncogenes (OGT, FoxM1, PI3K, c-Myc) and the presence of tumor suppressors (p53 and p21).

The graphs presented in Figs. 6–7 are sub-graphs, extracted from the complete state graph and show step-by-step changes towards progression of cancer and recovery, respectively. Since the complete state graph obtained from qualitative modeling is too complex to analyze each trajectory manually, we used the idea of betweenness centrality to identify important trajectories.

HyTech results

Hybrid modeling was carried out using the HyTech (HYbrid TECHnology) tool (Henzinger, Ho & Wong-Toi, 1997). The Bio-LHA of the significant cycle [(0, 1, 0, 0, 0, 0, 0, 1, 0) → (0, 1, 0, 0, 1, 0, 0, 1, 0) → (0, 1, 0, 0, 1, 1, 0, 1, 0) → (0, 1, 0, 0, 0, 1, 0, 1, 0) → (0, 1, 0, 0, 0, 0, 0, 1, 0)] defines invariants and clock rates for each qualitative state in the cycle (Fig. 8). The invariance kernel for the cycle (HyTech code in Supplemental Information 6) is composed of three conjuncted delay constraints, (Table 2). These delay constraints define necessary and sufficient conditions in a way such that the resulting trajectories maintain cyclic stability (homeostasis). If these delay constraints are violated, the trajectories will deviate from the significant cycle and may follow a path to the deadlock state. In Table 3, a pair-wise matrix represents an association between these delay constraints with the help of logical relations (≤, <, >, ≥, =).

The delay constraints from the start state (1,0,0,0,0,0,1,0,0) to the deadlock state (1,0,1,1,1,1,1,0,1) are computed (HyTech code in Supplemental Information 7) and presented in the form of a relation matrix in Table 4, which highlight important relations between the production and degradation rates of important proteins. The degradation delay of OGT is greater than the activation delays of MDM2, PI3K and p53, meaning that during progression towards deadlock state, depletion of OGT takes place at a much slower rate compared to MDM2, PI3K and p53 (Table 4). Similarly, production delay of c-Myc remains higher than PI3K and p53-MDM2 circuitry, thus reinforcing the results of our qualitative modeling which suggest that once c-Myc and OGT form a positive feedback loop, the overall system eventually moves to the deadlock state.

Qualitative dynamics

The qualitative model, rendered as a state transition graph (Fig. 5), highlights important behaviors in the form of trajectories to represent temporal evolution of the overall system from one qualitative state to another. These behaviors mainly include cycles that represent homeostasis, stable states (deadlock), and several bifurcation states from where the system can evolve either in the direction of tumor invasion or recovery. Three bifurcation states [(1,0,1,1,0,0,1,0,0), (1,0,1,1,1,0,1,0,0), and (1,0,1,1,1,1,1,0,0)] were observed to lead to both a deadlock state (1,0,1,1,1,1,1,0,1) and a typical reset state (0,0,0,0,0,0,0,0,0). Network analysis carried out using Cytoscape (Shannon et al., 2003) helped to identify important trajectories based on betweenness centrality and highlight step-by-step alterations that systematically lead to cancer metastasis from the starting state. Here, we discuss two cases that focus on important trajectories of the qualitative model involved in tumor progression and recovery.

Case 1: cancer initiation and progression

Under normal physiological conditions, proto-oncogenes play basic roles in signaling pathways that control cellular growth (Pall, 1981). Activation of a proto-oncogene into an oncogene, through gain-of-function mutations, increases the expression of these proteins and leads to alterations in signaling pathways, increases glycolytic flux through the HBP, and elicits an inflammatory response (Ma, Vocadlo & Vosseller, 2013; Karin, 2009; Fardini et al., 2013). Important trajectories that originate from a starting state (1,0,0,0,0,0,1,0,0) show increased expression of the pro-inflammatory NF-κB pathway and O-GlcNAc transferase (OGT) (Fig. 6). Supplemental Information 9 highlights changes in the expression level of genes along each transition that lead to the deadlock state. Enhanced OGT expression is considered an indicator of metabolic switch from oxidative phosphorylation to glycolysis, and can occur under oxidative stress. Cellular response to this stress has been shown to activate a p53 transcriptional response (Gambino et al., 2013), leading to a qualitative state (1,0,0,0,1,0,1,0,0). The subsequent trajectories, mapped in Fig. 6, show that p53 is over-expressed in several states that lead to a deadlock state, thereby reinforcing the growing evidence that suggests the divergent role of p53 in response to increased cellular metabolism (Puzio-Kuter, 2011; Maddocks & Vousden, 2011).

Despite its original classification as a tumor suppressor gene, recent evidence is accruing to reveal p53 also carries oncogenic properties (Soussi & Wiman, 2015). The majority of p53 germline and somatic alterations are missense mutations which synthesizes a stable mutant p53 protein that accumulates in the nucleus of tumor cells and can result in an oncogenic phenotype (Dittmer et al., 1993). Indeed, over-expression of p53 has been reported in various breast cancer studies, and induced a metabolic shift toward glycolysis (Won et al., 2012). Similar findings have also been reported in other tumor types, including cancers of the colon, cervix, and pancreas (Al-Khayal et al., 2016; Rajeshkumar et al., 2015; Hernández-Reséndiz et al., 2015; Kruiswijk, Labuschagne & Vousden, 2015).

In addition to accumulation of p53, significant increase in the expression of NF-κB has been reported in different tumor types, including breast and prostate cancers (Arora et al., 2014; Mak et al., 2015). Recent studies have shown the anti-apoptotic properties of hyper O-GlcNAcylation in tumor cells and the contribution of this post-translational modification for oncogenic activation of NF-κB in pancreatic cancers (Ma & Vosseller, 2013). In corroboration with these studies, our qualitative model demonstrates a sustained activation of NF-κB, which contributes to increased glycolytic flux and tumorigenesis.

The qualitative model reveals activation of FoxM1 and PI3K in response to continuous activation of OGT, thus leading to qualitative states (1,0,0,1,0,1,0,0) and (1,1,1,1,1,0,1,0,0). As expected, FoxM1 over-expression has been implicated in cancer. Although the mechanism of FoxM1 and OGlcNAcylation is poorly characterized, sentinel studies suggest that hyper O-GlcNAcylation of FoxM1 mediators in breast cancer prevent the degradation of FoxM1, to promote transformation of cells in breast cancer (Caldwell et al., 2010). From the initial state, betweenness centrality of states decreases as systems get closer to the deadlock state (Fig. 6). This drop in betweenness centrality is indicative of fewer chances for transition to recovery, particularly once FoxM1 and PI3K are over-expressed. Importantly, these results highlight the critical role of c-Myc to reach the deadlock state in the qualitative model. Prior to the activation of c-Myc, several bifurcation states exist for possible transitions to a typical reset or recovery state (0,0,0,0,0,0,0,0,0). However, qualitative modeling demonstrates that activation of c-Myc promotes the stability of hyper- O-GlcNAcylation. We observed that c-Myc activation forms a positive feedback loop with OGT, which plays a critical role in uncontrolled proliferation of tumor cells. Our findings are consistent with current literature that hypoxia-inducible factor (HIF), a transcription factor that activates aerobic glycolysis under cellular stress, cooperates with c-Myc to flip the metabolic switch and fuel glycolysis (Semenza, 2007; Dang et al., 2008; Kroemer & Pouyssegur, 2008). Together these proteins also upregulate glucose transporters, glycolytic intermediates, and induce angiogenesis in the tumor microenvironment to maintain glycolytic conditions (Kroemer & Pouyssegur, 2008).

Constitutive activation of c-Myc and p53 mutation contribute to uncontrolled cellular proliferation associated with upregulated glycolysis and metabolic re-programming in tumors. Further, activation of NF-κB downregulates oxidative phosphorylation in various tumor types (Markert, Levine & Vazquez, 2012). Together, our findings further characterize the critical roles of these oncogenes and tumor suppressor genes in support cancer progression through the regulation of biological networks.

Case 2: recovery from advanced tumor progression

In our qualitative model, important trajectories involved in cancer recovery suggest that the p53-Mdm2 circuit undergoes a series of cycles first to restore the system back to recovery state and subsequently, maintain homeostasis. These cycles are depicted as strongly connected components (SCC) in Fig. 7 (SCC1, SCC2, and SCC3). Supplemental Information 10 highlights change in expression level of genes along each transition that lead to the recovery state. In Fig. 7, the initial state (1,0,1,1,0,0,1,0,0) is a bifurcation state characterized by high expression of oncogenes and increased glycolytic flux. The initial state may lead to both recovery and deadlock states depending on genetic alterations regulating the expression of genes involved in signaling. The qualitative model shows that, in response to oncogene activation and increased glycolytic flux, p53 remains constitutively active in several successor states. Increased expression of p53 mediates the down-regulation of OGT by inhibiting c-Myc, thus leading to a state (1, 0, 1, 1, 0, 0, 0, 0, 0). Subsequently, down-regulation of PI3K is also triggered by p53- mediated inhibition through p21 and low expression of OGT, resulting in a qualitative state (1, 0, 1, 0, 1, 1, 0, 0, 0). In subsequent trajectories, the p53-Mdm2 circuit acts as a repair mechanism to systematically reduce the expression of oncogenes through a series of oscillations.

• The first cycle (SCC-1) shows an oscillation of the p53-Mdm2 circuit, while maintaining increased expression of FoxM1 throughout the cycle.

• In the second cycle (SCC-2), the p53-Mdm2 circuit oscillates to down-regulate the expression of FoxM1 until the system reaches to recovery state (0,0,0,0,0,0,0,0,0).

These findings are in agreement with previous experimental studies, which suggest that p53 levels showed a series of pulses in response to DNA damage. Uri Alon, in 2000, first reported that the p53-Mdm2 circuit show dampened oscillations in irradiated breast cancer cells (Bar-Or et al., 2000). Later studies confirmed these results by showing that the p53-Mdm2 circuit undergoes a series of pulses at regular intervals (Lahav et al., 2004; Lahav, 2009). The results presented in our study illustrate the role of the p53-Mdm2 circuit in a series of oscillations that lead to recovery state, consistent with previous studies (Poltz & Naumann, 2012; Abou-Jaoudé, Ouattara & Kaufman, 2009). A similar behavior is exhibited by our qualitative model in the form of two important cycles, shown as SCC-1 and SCC-2 (Fig. 7) before reaching a recovery state. However, this model does not provide any information about the number of iterations or time spent within each cycle. Depending on the extent of DNA damage, two scenarios are possible: (1) the recovery state (0,0,0,0,0,0,0,0,0) may also serve as an unperturbed stable state where the p53 level remains low; or (2) the system reaches a ‘limit cycle’/attractor (SCC-3) where it continues to oscillate indefinitely with constant time period and amplitude. In practice, the systems comprised of negative feedback loops, like circadian rhythms, are fully capable of producing sustained oscillations to maintain homeostasis.

Hybrid modeling

Cellular metabolism and intracellular signaling converge into a complex regulatory network that is regulated by key interactions. Importantly, these interactions that regulate these pathways vary in speed. While changes in gene expression occur at a slower rate, post-translational protein modifications tend to occur rapidly (Chubukov et al., 2014). In our study, the behavior of the p53-Mdm2 circuit is dependent on the time delay between p53-dependent induction and Mdm2-controlled repression. The delay constraints computed using HyTECH for the significant cycle [(0, 1, 0, 0, 0, 0, 0, 1, 0) → (0, 1, 0, 0, 1, 0, 0, 1, 0) → (0, 1, 0, 0, 1, 1, 0, 1, 0) → (0, 1, 0, 0, 0, 1, 0, 1, 0) → (0, 1, 0, 0, 0, 0, 0, 1, 0)] with lowest betweenness centrality, serves as an important attractor. It represents normal homeostasis characterized by low expression of oncogenes and p53-Mdm2 oscillations. Therefore, it is important to know the necessary and sufficient conditions that limit the system to maintain a homeostatic behavior.

The first state observed in the cycle (0,0,0,0,0,0,0,0,0) represents an initial configuration of the system. This state also represents a typical reset state after recovery. Different trajectories have been outlined that lead to this state (Fig. 7). The results of qualitative modeling (Fig. 7) show that once the system reaches the reset state, it enters into an attractor (SCC3). This cycle [(0, 1, 0, 0, 0, 0, 0, 1, 0) → (0, 1, 0, 0, 1, 0, 0, 1, 0) → (0, 1, 0, 0, 1, 1, 0, 1, 0) → (0, 1, 0, 0, 0, 1, 0, 1, 0) → (0, 1, 0, 0, 0, 0, 0, 1, 0)], characterized by p53-Mdm2 oscillatory behavior and low expression levels of OGT and oncogenic proteins, represents the normal homeostatic behavior of the overall system. Moreover, network analysis reveals that this cycle has the lowest betweenness centrality among all the cycles in the qualitative model, which makes it an attractive cycle. Once the system enters into this cycle (attractor), it tends to limit itself only within the cycle. Therefore, it is important to compute delay constraints that enforce the system to maintain this cyclic behavior.

Table 2 presents the invariance kernel of this important cycle. It is composed of four conjuncted delay constraints, which remain true within the cycle. Apparently, ($d_{\mathrm{p53}}^{-}\ge d_{Mdm2}^{+}$) is the most significant constraint which states that the degradation delay of p53 is greater than the production delay of Mdm2. That is, the rate of p53 synthesis must be greater than or equal to the production rate of Mdm2 for the system to maintain homeostasis. Experimental studies corroborate that a long time delay between the increase in p53 and the increase in Mdm2 would lead to oscillatory behavior (Tyson, 2004; Ciliberto, Novák & Tyson, 2005). Additional work demonstrated that the amplitude of oscillations in the p53-Mdm2 negative feedback loop are more variable than the period as a result of low-frequency noise in rates of protein production (Geva-Zatorsky et al., 2006). Fluctuations in the behavior of protein circuits produces biological response variations even between individual cells. In our study, we elucidated differences in the rate of cellular functions and identified the requirements needed for a system to maintain homeostasis. Taken together, these results provide a deeper understanding of the modulation of biological networks that play critical roles in tumorigenesis.