The mTOR signaling pathway (Figure 1) is initiated by insulin, insulin-like growth factor 1 (IGF1) and Ras (Laplante and Sabatini, 2012) along with others. The binding of insulin with insulin/IGF1 signaling (IIS) receptors causes its autophosphorylation with subsequent recruitment of insulin receptor substrate (IRS) with its cytosolic domain. Activated IRS then activates several downstream effector proteins including PI3K.

The role of PI3K pathway in metabolism and growth is well-established and its dysregulation could result in certain metabolic disorders and cancers (Carracedo and Pandolfi, 2008). In addition to the activation of JNK pathway (Vivanco et al., 2007) that down-regulates PTEN transcription and promotes cellular proliferation by hindering apoptosis, it is also involved in the activation of Akt (Carracedo and Pandolfi, 2008). PI3K phosphorylates and converts phosphatidylinositol (4,5)-bisphosphate (PIP2) to phosphatidylinositol (3,4,5)-trisphosphate to (PIP3) (Engelman et al., 2006; Manning and Cantley, 2007). This event is important for the activation of phosphatidylinositol dependent protein kinase 1 (PDPK1) that eventually stimulates Akt. The activation of Akt is achieved through phosphorylation at two sites i.e., Thr308 and Ser473. Activated PDPK1 phosphorylates Akt at position Thr308 whereas another protein mTORC2 phosphorylates it at Ser473 (Alessi et al., 1997; Sarbassov et al., 2005). Both these phosphorylation events are essential for the complete activation of Akt. Thobe et. al. examined the influence of PI3K pathway kinases on mTORC2 and found PI3K mediated regulation essential for mTORC2 recruitment and further activation Thobe et al. (2017).

Akt regulates several downstream proteins such as Tuberous Sclerosis proteins 1 and 2 (TSC1/TSC2) and FOXO. The heterodimer of TSC1 (hamartin) and TSC2 (tuberin) primarily inhibits the activity of mTORC1 via conversion of active Ras homolog enriched in brain (Rheb)-GAP into inactive Guanosine diphosphate (GDP)-bound Rheb (Inoki et al., 2003; Tee et al., 2003). Akt phoshorylates and inhibits TSC1/TSC2 in order to activate mTORC1 (Inoki et al., 2002; Manning et al., 2002; Potter et al., 2002; Roux et al., 2004; Ma et al., 2005). Akt also elevates the expression of mTORC1 indirectly through the inhibition of FOXO (Guertin et al., 2006; Chen et al., 2010; Zoncu et al., 2011). FOXO acts a homeostatic regulator of cellular energy production and consumption processes under energy stress conditions. Another role of FOXO is to increase the activation of Rictor (a major component of mTORC2) and subsequently mTORC1.

PI3K signaling is mainly buffered through PTEN (Carracedo et al., 2008). PTEN serves as a tumor suppressor and mostly found mutated in its phosphatase domain (Eng, 2003) in several cancers (Li and Sun, 1997; Steck et al., 1997) causing overactive PI3K signaling. PTEN hydrolyzes phosphatidylinositol (3,4,5)-trisphosphate (PIP3) to phosphatidylinositol (4,5)-bisphosphate (PIP2) (Figure 1). In this way, PTEN inhibits PIP3 dependent downstream signaling events like membrane recruitment and activation of AKT to prevent cell growth and proliferation. Hence, PTEN holds critical position in maintaining homeostasis through the inhibition of oncogenic transformation (Carracedo and Pandolfi, 2008).

Finally, the activated mTORC1 activates several downstream effectors mainly eukaryotic translation initiation factor 4E (eIF4E)-binding protein 1 (4E-BP1) and S6 kinase 1 (S6K1), that represses autophagy and promotes protein synthesis. S6K1 phosphorylates rictor causing mTORC2 disassembly (Dibble et al., 2009; Julien et al., 2010; Treins et al., 2010) and also degrades IRS proteins (Harrington et al., 2004; Shah et al., 2004) which dampens the PI3K mediated signaling cascade.

Mutations in the genes encoding for the proteins of mTOR associated BRN can lead to different types of cancers including sporadic cancers, hamartoma syndromes or phakomatoses, cowden syndrome (PTEN), neurofibromatosis (NF1, NF2) and peutz–Jeghers syndrome (LKB1) (Menon and Manning, 2008). Deregulation of entities in mTOR associated BRN result in certain other complications like insulin resistance and type 2 diabetes. Moreover, mTOR pathway can also be hyper-stimulated (e.g., in adipose tissues) under situation of excessive nutrients that can ultimately lead to the same complications (Um et al., 2004; Khamzina et al., 2005; Tremblay et al., 2007).

Gene expression is a complex process and its regulation determines the overall cellular dynamics (De Jong, 2002). Computational techniques in systems biology facilitate to explore the role of genes, proteins and overall dynamics of the system (Glass and Kauffman, 1973). Qualitative modeling framework is one of the established methods to analyze gene expression dynamics (Thomas, 1978; Thomas and d'Ari, 1990; De Jong et al., 2004) in the form of biological regulatory networks (BRNs) (Lewin, 2000). A BRN is modeled by a directed graph where vertices represent biological entities e.g., DNA, RNA, proteins and other biological molecules whereas edges correspond to regulatory interactions (i.e., activation and inhibition) (Bernot et al., 2007). The design of the study is illustrated in Figure 2.

The main objective of this study is to build a refined computational model of mTOR regulation that could predict therapeutic targets to inhibit the progression of cancer. A BRN of mTOR and its interacting proteins (PI3K, PTEN, mTORC2, Akt, mTORC1 and FOXO) has been abstracted from the pathway (Figure 1) in order to explore the dynamics based on the logical formalism of René Thomas (Peres and Jean-Paul, 2003; Bernot et al., 2004, 2007). The unknown parameters in the logical model are inferred based on biological observations formally encoded as CTL (Computational Tree Logic) formulas in SMBioNet (Selection of Models of Biological Networks) tool (Mcadams and Shapiro, 1995). The qualitative model (State Graph) of the BRN infers the dynamics such as homeostasis in the form of qualitative cycles and stable behavior in the form of stable state (SS). The most significant qualitative cycle is selected based on the centrality values of the qualitative states in the model. A linear hybrid automaton of the selected cycle is constructed using Hytech model checker (Henzinger et al., 1997) with new parameters for production and degradation time delays. Hytech inferred the values of these parameters in the form of linear constraints. These constraints are further analyzed to infer the pairwise relations between any possible pair of delays of genes. These relations show only one significant relation between AKT and mTORC2 in terms of delays. The analysis of the hybrid model provides a formal proof that during homeostasis the inhibition time delay of Akt is less than the inhibition time delay of mTORC2. This enforces that in homeostasis Akt is degraded with a faster rate than mTORC2 which suggests that the inhibition of Akt along with the activation of mTORC2 may be a better therapeutic strategy for the treatment of cancer.

The signaling pathway shown in Figure 1 is further reduced to a BRN shown in Figure 3 by the reduction rules described in (Naldi et al., 2009; Saadatpour et al., 2013). These rules have already been applied to reduce the TLR4 and JAK/STAT signaling pathways to a BRN with all possible regulatory feedback circuits (Paracha et al., 2014). The abstracted mTOR-associated BRN is composed of six proteins which are PI3K, PTEN, mTORC2, Akt, mTORC1, and FOXO.

Biological regulations (production and degradation) are subjected to expression levels of entities in BRN. An entity p1 activates or inhibits another entity p2, at a specific threshold. A qualitative threshold can be described as a discrete level (first, second, third etc.). René Thomas proposed a modeling framework which assumes qualitative thresholds and parameters to derive the dynamics of a BRN. Several methods are in use to model the behavior of biological systems (Peres and Jean-Paul, 2003; Bernot et al., 2004, 2007). Continuous modeling frameworks based on ordinary and partial differential equations are widely used. These frameworks rely on precise quantitative values, which in many cases are not known. This limitation led to the development of qualitative modeling framework. Kauffman et al., introduced a logical formalism based on Boolean logic where each entity was considered as “ON” (1) or “OFF” (0) to represent its activation or inhibition, respectively (Kauffman, 1969, 1993; Somogyi et al., 1997). This approach was further extended to kinetic logic formalism by Thomas to incorporate multi-valued (0,1,2,3,…) expression levels of entities. Formal methods such as model-checking approach can help to infer the parameters of complex systems (Bernot et al., 2007). BRNs are complex systems and their parameters can be inferred with such approaches.

This study is based on the kinetic logical formalism developed by René Thomas (Thieffry and Thomas, 1995) to model the biological regulatory network (BRN) of mTOR using GENOTECH tool Ahmad (2009) (available at https://github.com/DrJamilAhmad/GENOTECH/blob/master/GenoTechE.jar). An important feature of kinetic modeling is positive or negative feedback circuits. An entity favors the activation of another entity in the BRN through positive feedback and is necessary to generate multi-stationarity (stable states), whereas an entity favoring the inhibition of another entity through negative feedback is a necessary condition to generate oscillatory behavior (homeostasis) (Thomas, 1981). Number of studies performed on genetic networks that incorporated analysis of positive and negative feedbacks with formal methods can be found in Kauffman (1993), Somogyi and Sniegoski (1996), and Szallasi and Liang (1998). Formal definitions provided in Aslam et al. (2014) and Paracha et al. (2014) can be obtained for detailed description.

Qualitative dynamics of Thomas networks depend on the values of logical parameters which are unknown a priori. These parameters are used to render system dynamics as a directed state graph (discrete or qualitative model), which incorporates important behaviors such as cycles or stable states. The inference of biologically coherent parameters is an important aspect of qualitative modeling. In this direction, Bernot et al., introduced an approach to infer these logical parameters using a formal method approach called model checking. This approach is implemented in SMBioNet (Selection of Models of Biological Networks) tool. It performs an exhaustive enumeration of models and selects those set of parameters which are consistent with experimental observations expressed as temporal logic formulas. Similar parameter estimation approach (by using SMBioNet tool) has been employed to study qualitative behavior of several biological systems including immunity control mechanism in lambda phage network (Mcadams and Shapiro, 1995), pathogenesis and clearance mechanism of dengue virus (Aslam et al., 2014), MAL-Associated network of Cerebral Malaria (Ahmad et al., 2012) and the role of OGT in Cancer progression (Saeed et al., 2016).

Graph-theoretic approaches have been successfully applied on large protein networks (Barabasi and Oltvai, 2004; Stelniec-Klotz et al., 2012). The state graph can be further analyzed using network analysis techniques in terms of graph connectivity (Junker and Schreiber, 2011) by sorting it on the basis of maximum betweenness centrality (Freeman, 1977). The states with higher betweenness centralities represent higher chances of their occurrences. This may in terms of biological phenomenon represent the entities with frequent expressions. The qualitative states in the model with high betweenness centralities are compared to the rest of the state space in order to identify most favorable cycle (Tareen et al., 2015; Saeed et al., 2016).

René Thomas' framework provides useful insights into the discrete qualitative behavior of a biological system. However, naturally, the expression levels of proteins evolve in a continuous manner. Hybrid modeling combines discrete changes of a system with continuous changes (differential equation) in a single formalism Bio-Linear Hybrid Automaton (Bio-LHA) has been proposed for the hybrid modeling of qualitative BRNs (Ahmad et al., 2007). Bio-LHA uses time delays along with continuous variables (clocks) to compute production and degradation time of gene expressions. Production (δ⁺) or degradation (δ⁻) delay is the time required for a gene expression to reach from a lower level to a higher level or vice versa (Figure 4). In this approach, a clock variable (h) is associated with each entity which is initially set to zero and it evolves with rate 1 when the expression evolves. A clock is reset when it measures a production or degradation delay as shown in Figure 5. Hybrid model checking tool such as HyTech (Henzinger et al., 1997) can be used to infer the values of delays in the form of linear delay constraints for behaviors (paths toward stable states and cycles) observed in the qualitative model. Invariance kernel (Ahmad et al., 2007; Ahmad and Roux, 2010) represents cycles in the hybrid models which can also be characterized with delay constraints. These delay constraints are further converted into the relation matrix in order to find constraints between any two types of delays (production or degradation) of all entities. This modeling approach has been successfully applied to model a variety of BRNs (Ahmad et al., 2007, 2008, 2012; Ahmad, 2009; Aslam et al., 2014; Bibi et al., 2016).

To construct qualitative model of the mTOR associated BRN, SMBioNet tool has been employed to correctly estimate logical parameters according to biologicalobservations in literature (Chen et al., 2010; Zoncu et al., 2011; Laplante and Sabatini, 2012). This tool takes Computation Tree Logic (CTL) formulas representing biological observations and BRN as inputs and selects those sets of parameters which verify these formulas. In BRN modeling these parameters are used to incorporate behaviors in the form of paths, cycles and stable states (SS) as specified in the CTL formulas (Table 1). Formula 1 in Table 1 represents that in a particular homeostatic behavior (represented by CTL operator E) all entities with 0 expression levels after next qualitative state (represented by CTL operator X) finally reach the same expression levels in future (represented by CTL operator F). Formula 2 represents the biological observation that the overexpression of PTEN, FOXO, PI3K, mTORC1, and inactivation of Akt gradually leads to (represented by CTL operator ⇒) a steady carcinogenic state (represented by CTL operators F and G) with the overexpression of Akt, mTORC1, mTORC2, PI3K, and inactivation of PTEN and FOXO. The effect of PTEN inhibition on Akt/mTORC1 pathway eventually leads to a SS where FOXO and PTEN are not expressed. This behavior is encoded by Formula 3. On the basis of CTL formulas, SMBiogenerated eight sets of logical parameters for mTOR associated BRN (see Supplementary Files 1, 2).

The eight qualitative models of these sets were further analyzed for cycle(s) and SS(s) using GENOTECH tool. Almost all the parameter sets provide comparable results, revealing similar cycles and SS(s). First four models were selected on the basis of biological plausible SS (1, 0, 1, 1, 1, 0) representing the activation and inactivation states of entities in order of PI3K, PTEN, mTORC2, Akt, mTORC1 and FOXO, respectively. This SS represents the activation of PI3K, mTORC2, Akt and mTORC1 along with the inactivation of PTEN and FOXO. Subsequently, these 4 models were further compared for the logical parameter values which are coherent with biological observations. The set of logical parameters for this model (Supplementary File 3) is given in the last column of Table 2. The selected model (Figure 3) also shows eight cycles along with one SS (1, 0, 1, 1, 1, 0).

Several states lead the BRN directly into SS which represent critical divergence toward disease. All such states that eventually progress toward deadlock state do not possess functional PTEN or FOXO while having cellular proliferatory elements fully activated, e.g., (1,1,1,1,1,0), (1,0,1,1,1,1) as represented in (Figure 6). Thus, the down regulation of these tumor suppressor genes (PTEN and FOXO) bring this deadlock (SS) where no regulator is present to perform its function. The states (1,0,1,1,0,0), (0,0,1,1,1,0), (1,0,0,1,1,0), and (1,0,1,0,1,0) with temporary inhibition of PI3K, mTORC2, Akt, or mTORC1 are restrained in proceeding states resulting in their full and uncontrollable activation. This impact of tumor suppressor gene can be observed in cyclic state (1,1,1,1,1,0) that progresses into (0,1,1,1,1,0) (Figure 6) where PI3K is down-regulated through the inhibitory effect of PTEN to slow down further increase in cell mass and number. So the pre-occupation of PTEN is desired to recover the system into homeostatic state (0,1,1,1,1,0) that otherwise could divert into SS (1,0,1,1,1,0). The selected cyclic trajectory along with its respective constraints is given in Figure 8, specifies the stay conditions for each cyclic state and its violation would activate counter mechanism of autophagic inhibition by mTORC1 leading to cancer.

We applied exhaustive model-checking to validate the qualitative model of this study as proposed in Ben Abdallah et al. (2015). In this approach, the authors present a logical approach (using the Answer Set Programming language (ASP) Baral, 2003) to simulate and exhaustively analyze the dynamics of multivalued biological regulatory networks. The ASP method searches the attractor basins (stable states) which the region from which it is not possible to exit. By translating the model of Figure 3 to an automata network (Supplementary Files 6–9) and giving it as an input to the method of Ben Abdallah et al. (2015), we found that the set of all the attractor basins is reduced to a single stable state: (PI3K = 1, AKT = 1, PTEN = 0, MTORC2 = 1, FOXO = 0, MTORC1 = 1). This result is effectively coherent with the qualitative model given in Figure 6.

Since the model shows eight cycles therefore it is important to identify the most probable biological cycle. Thus, on the basis of betweenness centrality a cycle was computed by using Cytoscape tool (Shannon et al., 2003) that sorts all the states on the basis of their betweenness centralities (Freeman, 1977; Tareen et al., 2015), as presented in Figure 7 (Supplementary File 4). The nodes with larger diameter represent states with higher betweenness centrality. The cycle with maximum betweenness centrality: (1, 1, 1, 0, 0, 1) → (1, 1, 1, 1, 0, 1) → (1, 1, 1, 1, 0, 0) → (1, 1, 1, 1, 1, 0) → (0, 1, 1, 1, 1, 0) → (0, 1, 1, 0, 1, 0) → (0, 1, 0, 0, 1, 0) → (0, 1, 0, 0, 1, 1) → (0, 1, 0, 0, 0, 1) → (0, 1, 1, 0, 0, 1) → (1, 1, 1, 0, 0, 1)shows oscilation of all entities except PTEN. The cycle reveals that the constant activation of PTEN is required for homoeostasis. Any diversion from this cycle would either lead toward carcinogenic SS (1,0,1,1,1,0). All the cyclic trajectories represented in Figure 6 show expression of PTEN that positively regulates PI3K and enforces the model to remain in homeostasis. On the contrary, uncontrolled expression of PI3K and subsequent stimulation of cellular proliferative machinery mainly Akt and mTORC1 would either cause diabetic disorders or more severe circumstances like oncogenesis.

The cycle in Figure 7 shows homeostatic biological regulation of entities in the form of the switching of their low and high expression levels. In the cycle, the stable high expression of PTEN gene is revealed as a mandatory condition in all the states to maintain homeostasis (represented by 1 expression level in all cyclic states) while the expression of other entities oscillate in relation to each other. The Bio-LHA model in Figure 8 of this cycle was implemented in HyTech tool (Supplementary File 5) in order to predict its underlying causality relations of delays by analyzing the invariance kernel (Ahmad et al., 2008). The invariance kernel represents a set of viable cyclic trajectories in the state space of the Bio-LHA. Delay constraints characterizing the invariance kernel of the selected cycle were computed by HyTech tool (Table 3). In Table 3, the notation π is used to represent the sum of production and degradation delays as period (Ahmad, 2009). Conjunctions of all these constraints (1–9) constitute a necessary and sufficient condition for the existence of the invariance kernel and hence the qualitative cycle. Violation of any constraint would result in a null invariance kernel and eventually the qualitative cycle will no more exist. It is therefore sufficient that all the constraints should be valid (true) for the existence of the invariance kernel or qualitative cycle. For example, in Table 3, constraint 1 (δFOXO+≤|δmTORC1-|+|δAkt-|) shows that the production (activation) of FOXO is required before the degradation of Akt and mTORC1 and hence constitute a necessary condition for cycle (homeostasis). Again constraint 2 establishes another necessary condition for the existence of the cycle that explains that the degradation of FOXO should occur earlier than the production of mTORC1 and Akt. Similarly, other remaining constraints imposes necessary conditions for the existence of the cycle.

From the delay constraints in Table 3, a relation matrix is derived in Table 4 containing pairwise relations of delays of the entities FOXO, PI3K, mTORC1, mTORC2, and Akt. A relation between a pair of delays with both ≤ and ≥ reveals that the cycle is desensitized to the violation of such constraints. The table contains only one pairwise constraint between the degradation delays of Akt and mTORC2. It is important to note that these relationships are enforcing homeostasis represented by the qualitative cycle. From the results, it can be implicated that if cellular systems tends to escape homeostasis it may lead toward pathogenesis.

The only pairwise relation between δAkt- and δmTORC2- (Table 4) reveals a significant property of the selected homeostatic cycle where the degradation delay of Akt is less than or equal to the degradation delay of mTORC2. In other words, the degradatin of AKT occurs at faster rate than the degrdation of mTORC2. Interestingly, this is the only observed property of the cycle that enforces its existence dramatically. Thus, this constraint provides a governing rule that if violated may bifurcate the trajectory toward stable steady state (1,0,1,1,1,0) (represented by red colored state in Figure 6).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.