Journal of Theoretical Biology 355 (2014) 40–52

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Systematic modeling for the insulin signaling network mediated by IRS1 and IRS2 Can Huang a, Ming Wu b, Jun Du a, Di Liu a,n, Christina Chan b,c,d,nn a

Department of Mathematics, Michigan State University, East Lansing, MI 48824, United States Department of Computer Science and Engineering, Michigan State University, East Lansing, MI 48824, United States c Department of Chemical Engineering and Materials Science, Michigan State University, 2527 EB, East Lansing, MI 48824, United States d Department of Biochemistry and Molecular Biology, Michigan State University, East Lansing, MI 48824, United States b

H I G H L I G H T S







Construction of a systematic model for insulin signaling network mediated by IRS1 and IRS2. Current knowledge and mathematical models of insulin signaling pathways are integrated. Simulation shows how the wiring of the network determines different functions of IRS1 and IRS2. Sensitivity analysis identifies essential regulators for the signaling process.

art ic l e i nf o

a b s t r a c t

Article history: Received 16 May 2013 Received in revised form 26 February 2014 Accepted 19 March 2014 Available online 1 April 2014

The hepatic insulin signaling mediated by insulin receptor substrates IRS1 and IRS2 plays a central role in maintaining glucose homeostasis under different physiological conditions. Although functions of individual components in the signaling network have been extensively studied, our knowledge is still limited with regard to how the signals are integrated and coordinated in the complex network to render their functional roles. In this study, we construct systematic models for the insulin signaling network mediated by IRS1 and IRS2, through the integration of current knowledge in the literature into mathematical models of insulin signaling pathways. We hypothesize that the specificity of the IRS signaling mechanisms emerges from the wiring and kinetics of the entire network. A discrete dynamic model is first constructed to account for the numerous dynamic features in the system, i.e., complex feedback circuits, different regulatory time-scales and cross-talks between pathways. Our simulation shows that the wiring of the network determines different functions of IRS1 and IRS2. We further collate and reconstruct a kinetic model of the network as a system of ordinary differential equations to provide an informative model for predicting phenotypes. A sensitivity analysis is applied to identify essential regulators for the signaling process. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Discrete model Dynamic simulation Kinetic model

1. Introduction Liver is the central organ of glucose, lipid and energy metabolism in mammalian systems. A major function of the liver is to maintain glucose homeostasis under different physiological conditions. In the postprandial state (after food intake), liver cells

n Correspondence to: Department of Mathematics, Michigan State University, D217 Wells Hall, East Lansing, MI 48824, United States. Tel.: þ 1 517 353 8143; fax: þ1 517 432 1562. nn Corresponding author at: Department of Chemical Engineering and Materials Science, Michigan State University, 2527 EB, East Lansing, MI 48824, United States. Tel.: þ 1 517 432 4530; fax: þ 1 517 432 1105. E-mail addresses: [email protected] (D. Liu), [email protected] (C. Chan).

http://dx.doi.org/10.1016/j.jtbi.2014.03.030 0022-5193/& 2014 Elsevier Ltd. All rights reserved.

respond to elevated levels of serum insulin by promoting glycogen synthesis and inhibiting glycogenolysis and glyconeogenesis, while in the fasted state gluconeogenesis and other catabolic processes are activated to maintain serum glucose level and support energy consumption of the body, which is also mediated by insulin signaling (Gribble, 2005; Saltiel and Kahn, 2001). It has been shown that insulin signaling pathways play essential roles in such hepatic function to consume, store and produce glucose. Impaired insulin sensitivity (insulin resistance) or dysregulated insulin response in liver cells contributes to the development of type 2 diabetes (Michael et al., 2000). The hepatic insulin signaling is mediated by adapter proteins insulin receptor substrate (IRS), which can be phosphorylated on tyrosine residues by insulin receptors. IRS transmits the signal to downstream kinase cascades that regulate the expression of genes

C. Huang et al. / Journal of Theoretical Biology 355 (2014) 40–52

involved in the control of glucose and lipid metabolisms (White, 2003). There are two major forms of IRS: insulin receptor substrate-1 (IRS1) and insulin receptor substrate-2 (IRS2), both are abundantly expressed in normal liver cells (in human or mice) but down-regulated in diabetic patients and mice. Indeed, reduced IRS protein expression or reduced tyrosine phosphorylation of IRS in liver cells are hallmarks in the development of insulin resistance and type 2 diabetes (Gual et al., 2005; Zick, 2004). Essential downstream kinase cascades activated by IRS have been extensively studied, which involves the phophatidylinositol 3-kinase (PI3K)–AKT pathway (Brunet et al., 1999) and the phosphorylation of transcription factor forkhead box O1 (FoxO1) (Puigserver et al., 2003). Taniguchi et al. (2005) used short hairpin RNA to specifically knock down IRS1 and IRS2 in mouse liver, and showed that IRS1 and IRS2 work together in a complementary fashion to maintain PI3K and FoxO1 activity. Single knock-down of IRS1 or IRS2 did not cause significant diabetic phenotype in mice while a concomitant knock-down resulted in metabolic disorders, including hyperglycemia, hyperinsulinmia, insulin resistance and glucose intolerance (Taniguchi et al., 2005). Although IRS1 and IRS2 are highly homologous and mediate insulin signaling through the PI3K–AKT signaling pathway, they are suggested to assume distinctive functional roles in insulin signaling (Fritsche et al., 2008; Haeusler and Accili, 2008). This is evident by the phenotypic differences observed in knock-out mice (Kubota et al., 2008), whereby the liver specific IRS1 knock-out developed insulin resistance only after refeeding but not during fasting, while the IRS2 knock-out mice exhibited insulin resistance in the fasted state. These experiments suggest that insulin signaling adapts to different physiological functions through IRS1 and IRS2 signalings with different responses, but it is unclear how these homologous IRS proteins respond differently, given that they are mediating essentially the same PI3K–AKT pathway. There should be different features in the kinetics of their underlying network that drive the specificity, which is difficult to evaluate when only individual interactions or components of the network are monitored. Individual components in the insulin signaling network have been extensively studied in the last decade with their functional roles characterized through bio-chemical, molecular, genetic and biophysical approaches (Fritsche et al., 2008). Nevertheless our knowledge is still limited with regard to how the signals are integrated and coordinated in the complex network to adapt its function to different physiological conditions. Therefore, systematic modeling is required to provide a framework in which one can integrate the current knowledge of individual interactions and pathways to better understand the regulation of the insulin signaling network. In this study, we construct a systematic model for insulin signaling mediated by IRS1 and IRS2. This is different from previous models that mostly focused on a single pathway (e.g., MAPK; mitogen-activated protein kinase) (Hatakeyama et al., 2003), and IRS–PI3K–AKT (Sedaghat et al., 2002; Wu et al., 2009)). Instead, we integrate models of different pathways to demonstrate the difference in the functions of IRS1 and IRS2. We hypothesize that the specificity of the two IRS signalings emerges from the wiring and kinetics of the entire network. Based upon the current knowledge on this network, we reconstruct dynamic models by accounting for the following three important network features that could contribute to differential dynamics of IRS1 and IRS2 signaling pathways: (1) Complex feedback circuits within the network Phosphorylation of the tyrosine residues on IRS1 protein transmits insulin signaling, while serine/threonine phosphorylation of IRS1 inhibits the tyrosine phosphorylation, interrupts IRS1's interaction with insulin receptor (IR), or induces the

41

degradation of IRS1 (Paz et al., 1997; Pederson et al., 2001). Such inhibition of insulin signaling through serine/threonine phosphorylation of IRS1 can be catalyzed by IRS's downstream kinases, such as c-Jun N-terminal kinase (JNK) and I-κB kinase (IKK) (Zick, 2004, 2005), which have negative feedbacks. In contrast, we identified a positive feedback involving PKR (double-stranded RNA-dependent protein kinase), whose phosphorylation is down-regulated by insulin through the PI3K–AKT–PP1 pathway, while the phosphorylation of PKR inhibits insulin signaling by inducing the phosphorylation of IRS1 at serine 312 through kinases JNK and IKK (Wu et al., 2009; Yang et al., 2010). Insulin has also been shown to suppress the gene expression of IRS2 through a negative feedback (Hirashima et al., 2003; Zhang et al., 2001). The promoter region of the IRS2 gene contains an insulin response element, which can be recognized by the transcription factor FoxO1 (Zhang et al., 2001). High insulin concentration reduces the activity of FoxO1, thereby down-regulates IRS2 mRNA level in the liver (Zhang et al., 2001). As shown in our previous study in HepG2 cells, PKR also participates in this IRS–FoxO1 feedback by regulating the activity of FoxO1 (Yang et al., 2010). (2) Different time-scales of the regulations It has been shown that IRS1 is primarily regulated by posttranslational modifications such as tyrosine phosphorylation, serine phosphorylation and protein degradation, while IRS2 is regulated at the transcription level (Fritsche et al., 2008; Yang et al., 2010). in vivo studies in mice show that the mRNA level of IRS2 is increased during fasting and reduced after refeeding, while the mRNA level of IRS1 showed no significant change even after refeeding when IRS1 assumes a major role in response to insulin signaling (Kubota et al., 2008). Generally the regulations at the post-translation level occur quickly and are transient (limited by protein half-life) compared with regulations at the transcription level, which happens on a longer time-scale. Different regulatory time-scales of IRS1 and IRS2 are consistent with their physiological roles: a prompt response in the fed state controlled by IRS1 regulated at the post-translational level, and the long-term response during fasting controlled by IRS2 regulated at the transcriptional level. (3) Cross-talks between different pathways There are many downstream pathways in the insulin signaling network and the cross-talks between different pathways can play important roles in the regulation of the whole network. For example, our lab (C.Chan) has shown a potential cross-talk between Srchomology and collagen domain protein (Shc)– MAPK pathway and IRS–PI3K pathway through a JNK–PKR interaction. The phosphorylation level of PKR plays an important role in IRS–PI3K–AKT regulation through its feedback to IRS1, while PKR can also be regulated by the kinase JNK, which is downstream of the Shc–MAPK pathway. The Shc–MAPK pathway is another pathway that is activated in response to insulin; however it is not mediated by IRS. The activated insulin receptor interacts with protein Shc, which binds to proteins Grb2 and Son of Sevenless (SOS) to form a complex that activates the Ras-mediated MAPK pathway. The Shc–Grb2– SOS complex to activate Ras is modeled as “ShcGS” (Sasaoka and Kobayashi, 2000). Therefore, through the JNK–PKR interaction cross-talk occurs between two different insulin signaling pathways (Wu et al., 2009). Within the insulin signaling pathways mediated by IRS, the IRS–PI3K–AKT pathway represents a common route for both the IRS1 and IRS2 mediated signalings, even though different IRSs regulate different downstream target expressions. This suggests that cross-talk between the two signaling processes through the PI3K–AKT pathway may account for differential functions of IRS1 and IRS2 (Fritsche et al., 2008).

42

C. Huang et al. / Journal of Theoretical Biology 355 (2014) 40–52

By integrating current knowledge in the literature and models of insulin signaling, we reconstruct a network system consisting of these network features, and apply the discrete dynamic modeling approach developed by Wu et al. (2009) to study the dynamics of insulin signaling. We further collate and reconstruct a kinetic system for the network in the form of coupled ordinary differential equations to provide a more informative model for predicting phenotypes in the regulation of hepatic insulin signaling for better understandings of different functional roles of IRS1 and IRS2. We show the biological network of the discrete dynamic model in Section 2.1 and the kinetic model in Section 2.3. The simulation results and validation of the discrete model and the kinetic model are in Sections 2.2 and 2.3, respectively. The new biological discoveries we learned from the models are in Sections 2.4 and 2.5. All equations and parameters can be found in Section 3. We compare our models with previous models and discuss the advantages and limitations of different types of modeling in Section 4.

2. Results 2.1. Discrete dynamic model of IRS1 and IRS2 signalings Previously, we constructed a liver specific insulin signaling network mediated by IRS1 with essential feedbacks and cross-talks that were experimentally confirmed in HepG2 cells (Wu et al., 2009). We now extend the network to include IRS2 signaling. IRS2 can be activated (phosphorylation in tyrosine residues) by insulin receptors, to transmit the signal through the PI3K–AKT–PKR pathway. IRS2 is regulated at the expression level by transcription factor FoxO1. FoxO1 activity can be modulated by the phosphorylation of PKR through its regulation of the phosphatase PP1 (Yang et

al., 2010). The network is shown in Fig. 1 (details in IRS1 signaling network is shown in (Wu et al., 2009) and Fig. 1). A discrete modeling approach was applied to the IRS1 network to study the dynamics of the system (Wu et al., 2009). Here we introduce a novel feature into the model to account for the larger time-scale of transcriptional regulation of IRS2 as compared with kinase signaling (see Section 3).

2.2. Dynamic simulation shows different functional roles of IRS1 and IRS2 We use the discrete model to simulate dynamic responses of IRS1 and IRS2 signalings upon insulin stimulation under different physiological conditions. To mimic the physiological insulin serum level, we maintain a low insulin level at state 0 during fasting, and after refeeding the insulin signal jumps to the highest level (state 2), then gradually decreases along with time towards another fasted state. The activity of IRS1 and IRS2 upon such stimulation is shown in Fig. 2. IRS1 responds immediately after refeeding, while IRS2 takes a major role during fasting. The simulation shows that different IRS proteins assume the major role under different physiological states, which is controlled by the insulin level. Such a “dynamic relay” between IRS1 and IRS2 is consistent with in vivo experiments (Kubota et al., 2008). The discrete model is based on the network wiring or connections and assumes that the network topology/architecture defines the major dynamic characteristics of the system, thus our current simulation confirms that different dynamics and time-scales of regulations of IRS1 and IRS2 emerge from the regulatory logics in the complex network. However, since a discrete model is based solely on the qualitative information of the interactions in the network, it provides a rough approximation of the system by assuming similar kinetic

Fig. 1. A schematic representation of the IRS1 and IRS2 signaling network. IR is activated in response to insulin stimulation. IRS1 and IRS2 can be activated by IR through the phosphorylation of tyrosine residues on IRS proteins. Both IRS1 and IRS2 transmit insulin signal through PI3K pathways, which affect downstream PKR phosphorylation. PKR phosphorylation could induce the activation of FoxO1, which regulates the gene expression of IRS2. A detailed network of IRS1 signaling is also shown by Wu et al. (2009). “IRST” and “IRSS” in the figure represent the tyrosine phosphorylation and serine phosphorylation of IRS1, respectively. Colors of the edges describe different regulatory relationships. Arrow ends represent activation while blunt ends represent inhibition.

C. Huang et al. / Journal of Theoretical Biology 355 (2014) 40–52

43

Fig. 2. Simulation of IRS1 and IRS2 activity in response to changes in insulin signal. The simulation is based on the discrete dynamic model; the states of the components are averaged with 5000 different runs each with noises in the signaling process. Insulin level is set to Low (state 0) in the initial fasted state. During feeding, the insulin concentration is elevated immediately to a high (state 2) level and gradually decreases along with time, and finally reaches another fasted state (state less than 1). The plot shows the activity change of IRS1 and IRS2 in response to insulin level under different physiological states. The X-axis shows the “time steps” in the simulation of the discrete model.

constants for all reactions, which may lose subtle dynamics as well as introduce artifacts in the discrete model (Wu et al., 2009). It is possible that certain kinetic parameters are essential for the dynamics of IRS signaling, which would not be easily captured and evaluated by our discrete modeling approach. Another limitation of the discrete dynamic model is that there are only “time steps” as the unit of time in the model, which does not represent the actual time in the experiments. Therefore, we further collate mathematical models that describe different aspects of insulin signaling in the literature, and integrate them to reconstruct a model that combines current knowledge in the literature and accounts for various feedbacks, time-scales and cross-talks between different pathways. 2.3. Kinetic modeling of IRS2 and IRS2 signalings in liver cells We adopt and integrate the kinetic model of Shc–MAPK pathway by Fritsche et al. (2008), and the IRS–PI3K–AKT model by Sedaghat et al. (2002). The original model for IRS–PI3K–AKT pathway consists of regulations in different layers of insulin signaling, including receptor internalization, protein degradation and feedbacks to regulate IRS activity (Sedaghat et al., 2002), but it is not specific to the liver. The concentrations of the molecules in the model are not consistent with physiological conditions, e.g., a concentration of 10  13 M (or 0.0001 nM) for insulin receptor suggests that each cell has no more than a couple of the protein molecules, which is not likely the case in the liver. We address this shortcoming by modifying the feedbacks with the PKR mediated signaling pathways specific to liver cells that had been uncovered in our previous study (Wu et al., 2009; Yang et al., 2010). We connect this submodel into the kinetics of Shc–MAPK pathway by Hatakeyama et al. (2003), by introducing cross-talks through JNK (Wu et al., 2009; Yang et al., 2010). The transcriptional regulation of FoxO1 of IRS2 is also modeled with kinetic equations (see Section 3). Fig. 3 shows the major pathways that have been integrated in our kinetic model. We simulate our integrated ordinary differential equation model for the insulin signaling network, and compare the results with our experimental studies of short-term IRS activity and PKR phosphorylation level under insulin treatment (Wu et al., 2009;

Yang et al., 2010) (Fig. 5, compared with Fig. 3B in Wu et al. (2009)) which shows the time-series of IRS1 tyrosine and serine phosphorylations, and PKR activation upon insulin treatment in HepG2 cells.), thereby determine the missing parameters including the reaction rates for the feedbacks through PKR, JNK and FoxO1 regulations. The full list of parameters and their references are provided in Table 1. The equations and parameters define a kinetic model for the IRS1 and IRS2 signalings in liver. We also simulated the knock-out/inhibition of JNK (Fig. 6). In both models the inhibition of JNK enhances IRS1 activity (the tyrosine phosphorylation of IRS1), which is consistent with our experimental results (Wu et al., 2009; Yang et al., 2010; in particular, Fig. 5 in Wu et al. (2009)). The kinetic model is then used to simulate the dynamic responses of IRS1 and IRS2 signalings upon insulin stimulation under different physiological conditions. Fig. 4 shows that IRS2 activity (mRNA level) is elevated in the fasted state, while there is a transition between IRS1 and IRS2 activities upon refeeding. IRS1 is activated immediately after refeeding and the high level of insulin concentration inhibits IRS2. When the insulin signal gradually decreases and the system returns to the fasted state, IRS2 activates to assume a major role during fasting. The kinetic model confirms results of the network dynamics uncovered by the discrete model. In the discrete dynamic model, the downstream signal molecules are “activated” by upstream signaling molecules, without kinetic mechanisms of such activation. In contrast, the kinetic model provides such details. For example, the activation of ShGS complex by the insulin receptor involves the endocytosis of the receptor upon phosphorylation by insulin, the binding of Shc factor by the cytosolic fraction of the receptor, and the recruitment of Grb/SOS factors into an active complex. Each process here involves certain equilibrium of molecules between their different forms, e.g., the membrane-bound receptor and cytosolic receptor, the inactive Shc molecules and the receptor-bound Shc molecules, the Grb/SOS (GS) inactive form and the Shc/Grb/SOS (ShGS) active form. Similarly, the activation of Ras involves the GTP binding and hydrolysis process, as indicated in Table 1 and Fig. 3 (RasGTP, RasGDP), and the activation of PI3K involves the transition between different phosphorylation forms of phosphatidylinositol. All these detailed mechanism cannot be captured by the discrete

44

C. Huang et al. / Journal of Theoretical Biology 355 (2014) 40–52

Fig. 4. Simulations of IRS1 and IRS2 activity upon an insulin signal. The simulation is based on the kinetic model. Ordinary differential equations are simulated using the forward Euler method with step size 5  10  4 in MATLAB. The plot shows the activity change of IRS1 and IRS2 in response to insulin level under different physiological states.

shows the kinetic constant of MEK–JNK activation are the most sensitive in the system (Table 3), suggesting the cross-talk through MEK–JNK activation is the most important factor in the system. Thus the model predicts that inhibition of the JNK activation could effectively alter insulin response. Indeed, besides directly perturbing insulin substrates (e.g., the inhibition of the serine phosphorylation of IRS1; Sykiotis and Papavassiliou, 2001), down-regulation of JNK activity has been suggested to be an effective therapeutic target for the treatment of insulin resistance (Bennett et al., 2003). Therefore, our model could be used to identify important genes in the pathway that could be potential targets. 2.5. Characteristics of insulin signaling network learned from modeling and simulations

Fig. 3. Schematic representation of the processes modeled in the kinetic system. (A) IRS1 and IRS2 mediated insulin signaling pathways. The model includes the insulinreceptor binding, internalization and degradation of receptors, the phosphorylation of downstream kinase cascades through PI3K and PKR, and the transcriptional regulation of IRS2. Since many proteins involved in the process have different phosphorylation forms, the solid lines represent changes between different forms of a protein/protein complex, while the dashed lines represent regulatory relationships. We use arrow ends for activation and blunt ends for inhibition. (B) The MAPK–JNK pathway is activated independent of IRS adapters. Shc is activated by insulin receptors and initiates downstream kinase cascades upon insulin stimulation.

model, which describes only the ON/OFF (activated/inhibited) states of the components. Fig. 5 shows examples of the dynamic profiles that can only be captured by the kinetic model, namely the rapid phosphorylation of phosphatidylinositol upon insulin stimulation (Fig. 5D), and the decrease of RasGDP after insulin peaks (Fig. 5E). 2.4. Sensitivity analysis to identify key players in IRS1/2 signaling One advantage of the kinetic model is being able to identify pivotal components and key reactions by sensitivity analysis. We performed sensitivity analysis on the system with respect to all parameters and initial conditions (see Section 3). The analysis

The dynamic model integrates current knowledge of the IRS1/2 signaling network. Simulations of the dynamic model enhanced our understandings of how the network functions as a comprehensive system, which is beyond the ability of single experimental studies in which one or a few responses are evaluated upon perturbation to infer the underlying mechanism. Through the process of collating the information, modeling and simulation, we learned that the dynamics of the signaling network emerges from its network features: 2.5.1. Identification of the most important pathways that regulate insulin response Our model suggests the importance of feedback circuits through the serine/threonine phosphorylation of IRSs. The discrete dynamic model shows that the IRS responses change significantly when the negative feedbacks through IKK/JNK or the positive feedback through PRK is perturbed (Wu et al., 2009), however it is not able to distinguish the effect of every pathway given its allor-none assumption in signaling. In contrast, the kinetic model provides more details. In the sensitivity analysis on the kinetic constants (Table 3), JNK and Ras show the highest coefficient, suggesting that the activation/inactivation of JNK and Ras are the most sensitive in the system. In the sensitivity analysis on the initial conditions (Table 4), in addition to the insulin receptors (IR) and substrates (IRS) that were also shown in the discrete dynamic model, concentrations of PKR and Ras are the most significant

C. Huang et al. / Journal of Theoretical Biology 355 (2014) 40–52

45

Fig. 5. Simulations of IRS1, IRS2 and kinase activity under insulin treatment. This figure is obtained using the parameters listed in Table 1. (A) PKR is inhibited under insulin stimulation, and then gradually recovers. (B) Dynamics of total concentration of IRS2 in the fed state (high level of insulin) is shown. (C) This simulation is compared with our experimental results shown by Wu et al. (2009), Yang et al. (2010), specifically Fig. 3B in Wu et al. (2009) that shows the time-series of IRS1 tyrosine and serine phosphorylation, and PKR activation upon insulin treatment in HepG2 cells. The kinetic modeling results correspond to the experimental results reported by Wu et al. (2009), Yang et al. (2010) are shown in (D), (E), (F) and (G) through dynamics of PI(3,4,5)P3 and PI(3,4)P2, RasGDP, activated Akt and activated JNK, respectively.

46

C. Huang et al. / Journal of Theoretical Biology 355 (2014) 40–52

Fig. 6. Simulations of IRS1 and IRS2 activity under insulin treatment when JNK is knocked-out. (A) Responses of IRS1 and IRS2 to different levels of insulin signal with JNK knocked-out. (B) A comparison of tyrosine phosphorylated IRS1 before and after JNK knock-out.

coefficients. These results highlight the importance of the components PKR, JNK, and Ras in insulin signaling, among all the different feedback systems. In the wiring of the network, the components PKR, JNK and Ras are all local hubs coordinating multiple pathways and feedbacks. JNK integrates IRS-dependent and IRS-independent (Ras–MAPK) signalings. PKR receives signals from JNK and AKT pathways and mediates feedbacks through PP2A and IKK. Ras–Rac is the integration step of the ShGS pathway and PI3K signals. Therefore, points of cross-talks and feedbacks between signaling pathways could be potential therapeutic targets.

2.5.2. Different time-scales of signaling and transcription are the design principle underlying the IRS1/IRS2 switch under different physiological conditions The models that consider only the signaling process (involving receptor binding and kinase cascades, e.g., discrete model in Wu

Fig. 7. The effect of “delay” on the switch of IRS1 and IRS2. (A) Inhibit the transcription of IRS2 mediated by FoxO1. The threshold of [Active FoxO1] is increased from 2.5  10  7 M to 2.5  10  3 M. (B) Enhance the transcription of IRS2 mediated by FoxO1. The threshold of [Active FoxO1] is decreased from 2.5  10  7 M to 2.5  10  10 M.

et al. (2009), kinetic models in Hatakeyama et al. (2003) and Sedaghat et al. (2002)) are not able to show differences in functions of IRS1 and IRS2. Our discrete dynamic model is based on the model of Wu et al. (2009), but accounted for the transcriptional regulation of IRS2 with a delayed step. With such a small change we can then observe a switch between IRS1 and IRS2 in fasting and after refeeding (Fig. 2A). Only when we distinguish the time-scale of transcriptional regulation from other signaling processes mediated solely by protein–protein interactions could our discrete model explain the switching between IRS1 and IRS2. Thus we hypothesize that different time-scales in signaling and transcription are the key to understanding the switch. We tested the hypothesis in the kinetic model. Fig. 7 shows the altered switching behaviors upon change of the parameter in the FoxO1 mediated

C. Huang et al. / Journal of Theoretical Biology 355 (2014) 40–52

transcription. The delay of transcription eliminates the switching behavior upon insulin stimulation (Fig. 7A, 15–65 min). Interestingly, reducing the time-scale in transcription (enhance transcription) also results in a delay of the switch (Fig. 7B, the switch occurs at  40 min, later than in the normal state shown in Fig. 4 where it occurs at 20 min). Therefore, the efficiency of the switching behavior in the IRS response relies on the coordination of timescales in the transcription and signaling process. Overall, the modeling and simulation studies support our hypothesis that the dynamics and the specificity of the two IRS signalings emerge from the network features: the cross-talks and feedbacks in signaling, and the different time-scales between signaling and transcriptional regulation.

3. Method 3.1. Discrete dynamic modeling of IRS1/IRS2 signaling network The basic assumptions of the dynamic model are: (a) the wiring of the network determines the major dynamic characteristics; (b) the activity level of the components are represented by discrete status (ON or OFF); (c) given that modularity is a design principle of biological systems, isolated or individual modules of the system should be able to provide insight into the mechanisms that govern the system behavior; and (d) the cell variation is represented by fluctuations in their initial protein activity levels. The signaling network is formulated by a directed graph in which each directed edge represents a positive (activation) or negative (inhibition) regulatory interaction. We associate each component in the network with a three-state variable to represent its protein activity: 0 (lower than control, suppressed), 1 (control, normal level), and 2 (higher than control, activated). Time is modeled by regular intervals called time steps. Model simulation is essentially a series of state transitions that represent the changes in the states of each variable along the time steps. The transition rules are defined based on the regulatory interactions: if a component/protein A is positively regulated by B, when B is in the activated state (state 2), the state of A will shift up in the next time step (0-1, or 1-2), while for inhibitors the effect is opposite to that. The state of the system (states of individual components) is renewed based on the transition rule at each time step. The transition rules are based on activation or inhibition of a certain kinase cascade. Therefore, the state of the system reflects how the kinases are activated/inhibited upon insulin treatment and how long (i.e., how many intermediate steps) the cascade responds. The simulation averages multiple runs, where each starts from random initiation states and evolves asynchronously to represent possible heterogeneities in the cell population. Each simulation is independent with random initiation states of all components, except that the insulin level is turned ON. Averaging over 5000 simulations, the result shows the average evolving (temporal) profile of each component in the network in response to insulin signaling. Further details of the modeling approach can be found in Wu et al. (2009). Transcriptional regulations of FoxO1 on IRS2 expression are expected to respond much slower than other signaling processes on the post-translational level, thus we introduce a delay factor on the interaction such that the current level of IRS2 depends on a further previous state of FoxO1 activity, thereby the new simulation approach changes the original assumption of single-step Markov process Si ¼ f (Si 1) to a history-dependent state transition Si ¼f(Si 1, Si 2, … S1). In the discrete model, the state transition is in units of time steps and all reactions take one time step. In order to implement a delay in some reactions, intermediates are added between activator/inhibitor and their substrates, so that the state change of the substrate takes multiple time steps in response to the regulator. In practice, we

47

implement the idea by simply introducing multiple intermediates between FoxO1 and IRS2 to delay the signaling transduction. The dynamic model was implemented by customized MATLAB codes.

3.2. Kinetic modeling of IRS1/IRS2 signaling network The complete biological network subject to modeling is illustrated in Fig. 3. Different from the discrete dynamic model, we do not assume discrete status on the variables. More importantly, the dynamics are determined not only by the wiring of the network but also the kinetic processes/parameters in the network, which allows one to study more detailed mechanisms. We adapted some reaction from Hatakeyama et al. (2003) and Sedaghat et al. (2002) where the equations are the same except for several of the parameters. Hence, in this section, we focus on the newly introduced reactions that are described in some detail. For the complete set of reactions in the model, see Tables 1 and 2. (1) In the PKR feedback, PKR activates serine phosphorylation of IRS1. We described such activation as the regulatory effect of PKR upon the IRS serine phosphorylation represented by the following reactions: 0

k7

IRS1 -IRSSerP ; 1

0

k7

IRS1 ’ IRSSerP : 1

The PKR regulation is incorporated by adding to reaction rate 0 k  7 a multiplicative term 1 þ ðK max ½PKR=K M þ ½PKRÞ, which is equal to 1.0 in the case when PKR is not activated, and takes a maximal activation as PKR is 100% activated. Parameters Kmax and KM describe effects of maximum and half activations of PKR, respectively. When considering the AKT feedback through IRS serine phosphorylation, another multiplicative term 1 þ ðK max ½AKT=K m þ ½AKTÞ is further added. Therefore the differential equation for the reactions reads as 8 d½IRS1  > 0 0 > > ¼  k7 ½IRS1  þ k  7 ½IRSSerP ; < dt 1 SerP > > d½IRS1  ¼ k0 ½IRS   k0 ½IRSSerP ; > : 1 1 7 7 dt where 0

k7 ¼ k 

   ½IRp  þ ½IR2p  b1 ½PKRp  b2 ½AKT p  1 þ  1þ IRp 50 þ ½PKRp  50 þ ½AKT p 

IRp is the concentration of phosphorylated surface receptors after maximum insulin stimulation and k is the normal reaction rate assuming about 2 min half time of the dissociation. [PKRp] and [AKTp] are normalized to the total amount of PKR and AKT proteins, respectively. In the following, we will choose parameters such that b1 ¼ 50; b2 ¼ 2; IRp ¼ 8:97  0 10  9 and k  7 ¼ 0:516 (2) The phosphorylation of IRS2 is similar to that of IRS1, without the serine site inhibition. The reactions k7

IRS2 -IRStyrP 2 ;

″

k7

IRS2 ’ IRStyrP 2

″

where k  7 ¼ 1:4 lead to equation 8 d½IRS2  > > > ¼  k7 ½IRS2  þ k  7″ ½IRStyrP < dt 2 ; tyrP > tyrP > d½IRS2  ¼ k ½IRS   k > : 7  7″ ½IRS2 : 2 dt (3) As in Sedaghat et al. (2002), we assume the association of phosphorylated IRS2 with PI3K to occur with a stoichiometry of 1:1 and with the same reaction rates as that of IRS1/PI3K

48

Table 1 All species and parameters of IRS–PKR and MAPK signaling in our model. Process

Reactions

Rate equations

Parameters

Reference

Insulin-receptor binding (Eqs. (1)–(9))

Insulin þ Receptor2IR

k1 ½Insulin½Receptor  k  1 ½IR

k1 ¼ 6  107 ; k  1 ¼ 0:2

Sedaghat et al. (2002)

Insulin þ IRp 2I2Rp

k2 ½Insulin½IRp   k  2 ½I2Rp 

k2 ¼ k1 ; k  2 ¼ 102 k  1

Sedaghat et al. (2002)

IR2IRp

k3 ½IR  k  3 ½IRp 

k3 ¼ 2:5  103 ; k  3 ¼ k  1

0

0

Sedaghat et al. (2002)

0 0 0 k4 ½IRp   k  4 ½EN  IRp  k4 ¼ 2:1  10  3 ; k  4 ¼ k4 =10 0 0 k4 ½IRp   k  4 ½EN  IRp  k4 ½Receptor  k  4 ½EN_Receptor k4 ¼ 3:33  10  4 ; k  4 ¼ 9k4 2EN_Receptor (synthesis and k5  k  5 ½EN_Receptor k5 ¼ 1:67  10  18 ; k  5 ¼ 10k5 degradation) k6 ¼ 0:461 k6 ½EN_IRp  EN_IRp - EN_Receptor

IRp 2EN  IRp I2Rp 2EN  I2Rp Receptor2EN_Receptor

k6 ½EN_I2Rp 

IRS1 2IRSTyrP 1

k7 ½IRS1   k  7 ½IRSTyrP  1

IRS2 2IRSTyrP 2

00 k7 ½IRS2   k  7 ½IRSTyrP  2

IRS1 serine phosphorylation (Feedback: PKR, AKT, JNK) (Eq. (12))

IRS1 2IRSSerP 1

PI3K activation (Eqs. (13)–(16))

þ PI3K2 PI3K_IRSP1 IRSTyrP 2

IRS1=2 tyrosine phosphorylation (Eqs. (10) and (11))

Sedaghat et al. (2002) Sedaghat et al. (2002) Sedaghat et al. (2002)

″ k7 ¼ 4:16ð½IRp  þ ½I2Rp Þ=IRp; k  7 ¼ 1:4½PTP þ 109 a½IRSSerP ; k  7 ¼ 1:4; 1

Sedaghat et al. (2002), Guo et al. (2006)a

k7 ½IRS1   k  7 ½IRSSerP  1

½PTP ¼ 1  0:25½AKT p ; a ¼ 0:2; IRp ¼ 8:97  10  10    P P P P 0 Þ 0   k7 ¼ 2ð½IR  þ ½I2R 1 þ 50b1þ½PKR 1þ 50b2þ½AKT ; k  7 ¼ 0:516; b1 ¼ 50; b2 ¼ 2 8:9710  9 ½PKRP  ½AKT P 

Cedersund et al. (2008), Sedaghat et al. (2002)b,c

k8 ½IRSTyrP ½PI3K 1

k8 ¼ 7:0  1011 ; k  8 ¼ 10

Sedaghat et al. (2002)

0

0

 k  8 ½PI3K_IRSP1  IRSTyrP þ PI3K2 PI3K_IRSP2 2

k8 ½IRSTyrP ½PI3K 2

PI45P22PI345P3

 k  8 ½PI3K_IRSP2  k9 ½PI45P2  k  9 ½PI345P3

k9 ¼ ðk9s  k9b Þ

Guo et al. (2006)

½PI3K_IRSP1  ½PI3K total  þ k9b ;

PI34P22PI345P3

k10 ½PI34P2  k  10 ½PI345P3

k10 ¼ 3:08; k  10 ¼ 2:8

AKT activation (Eq. (17))

AKT2AKT P

k11 ½AKT  k  11 ½AKT P 

k11 ¼ 0:107ð½PI345P3  0:31Þ k  11 ¼ 3

PKR activation (Eq. (18))

PKR2PKR

P

ShGS complex formation (Eqs. (19)–(24)) Rac activation (via ShGS or PI3K)

IRP ðI2RP Þ þ Shc2IRSh IRSh2IRSh

P

P

P

k15 ¼ 5:4  104 ; k  15 ¼ 2:574

IRShGS2ShGS þ IRP

k16 ½IRShGS  k  16 ½ShGS½IRP 

k16 ¼ 7:2; k  16 ¼ 1:44  103

P

P

k17 ½ShGS  k  17 ½GS½Sh  P

P

Kholodenko et al. (1999) Kholodenko et al. (1999)b Kholodenko et al. (1999)b Kholodenko et al. (1999)b

k17 ¼ 6; k  17 ¼ 1:26  105 7

7

Hatakeyama et al. (2003), Kholodenko et al. (1999)b

V 18 ½Sh =ðK 18 þ ½Sh Þ

V 18 ¼ 1:02  10

RasGDP- RasGTP½ShGS

ke1 ½ShGS½RasGDP

ke1 ¼ 13:32; K e1 ¼ 1:81  10  10

Hatakeyama et al. (2003)b

RasGDP- RasGTP½PI3K

=ðK e1 þ ½RasGDPÞ kj1 ½PI3K½RasGDP

kj1 ¼ 13:32; K j1 ¼ 1:81  10  10

Hatakeyama et al. (2003)b,e

RasGTP -RasGDP

=ðK j1 þ ½RasGDPÞ V j2 ½RasGTP=ðK j2 þ ½RasGTPÞ

V j2 ¼ 1:73  10  8 ; K j2 ¼ 5:71  10  11

Hatakeyama et al. (2003)b

kj3 ½MEK½RasGTP =ðK j3 þ ½MEKÞ

kj3 ¼ 210; K j3 ¼ 3:17  10  7

Hatakeyama et al. (2003), Schoeberl et al. (2002)b Schoeberl et al. (2002)

P

MEK -MEK JNK-JNK

P

JNK P -JNK

; K 18 ¼ 3:40  10

Sedaghat et al. (2002)b,d

Sh 2Shc

MEKK activation (Eqs. (28) and (29)) MEK-MEK P

JNK activation cross-talk: PKR (Eqs. (30) and (31))

P

3

Kholodenko et al. (1999)

k14 ¼ 1:2  103 ; k  14 ¼ 3:0  102

k15 ½IRSh ½GS  k  15 ½IRShGS

P

Ras activation (Eqs. (25)–(27))

k13 ð½IRP  þ ½I2RP Þ½Shc  k  13 ½IRSh k13 ¼ 5:4  105 ; k  13 ¼ 36

IRSh þ GS2IRShGS ShGS2GS þ Sh

Sedaghat et al. (2002) P

k12 ¼ k12c ð1þ 10½JNK Þ k  12 ¼ k12c ð0:1 þ 2:99½AKT Þ k12c ¼ 6:9  10

k14 ½IRSh  k  14 ½IRSh 

P

Sedaghat et al. (2002)

P

k12 ½PKR  k  12 ½PKR 

Sedaghat et al. (2002)

k9s ¼ 1:39; k9b ¼ :125; k  9 ¼ 41:7;

kj4 ½MEK =ðK j4 þ ½MEK Þ

kj4 ¼ 3:48  10  9 ; K j4 ¼ 2:20  10  6

kj5 ½MEK½JNK=ðK j5 þ ½JNKÞ

kj5 ¼ 5:7  102 ; K j5 ¼ 1:46  10  4

kj6 ½JNK P =ðK j6 þ ½JNK P Þ

kj6 ¼ 1:8  10  8 ; K j6 ¼ 1:60  10  7

P

P

Hatakeyama Schoeberl et Hatakeyama Schoeberl et

et al. (2003), al. (2002)b,f et al. (2003), al. (2002)b

C. Huang et al. / Journal of Theoretical Biology 355 (2014) 40–52

EN_I2Rp - EN_Receptor

Sedaghat et al. (2002)

No direct or indirect evidence available thus the value is derived from model fitting. No evidence or estimations currently available, the values are borrowed from a different system and the original paper is referenced. The value of the rate constant is borrowed from PKC (Sedaghat et al., 2002), and we use the equation from Cedersund et al. (2008) to model the combinatory effect of AKT and PKR. We used the same form of equation for the effect of JNK, which is currently unknown. d The values of the rate constants of PKR activation are borrowed from PKC (Sedaghat et al., 2002). e The equations and parameters for Rac are borrowed from Raf (Hatakeyama et al., 2003), which is a similar GTP-binding protein that can activate MEKK. f MEKK activates both ERK and JNK, so the equations and parameters for JNK are borrowed from ERK (Hatakeyama et al., 2003). The same model can be used for JNK activation through PKR.

49

activations k8

IRStyrP þPI3K-PI3K_IRSP2 ; 2 k8

IRStyrP þPI3K ’ PI3K_IRSP2 : 2 The corresponding equation for involved reacting species with the newly introduced reactions is given by 8 d½IRStyrP > P 2  > > ¼  k8 ½IRStyrP > 2 ½PI3K þk  8 ½PI3K_IRS2 ; > dt > > < P d½PI3K_IRS2  P ¼ k8 ½IRStyrP 2 ½PI3K k  8 ½PI3K_IRS2 ; > dt > > > > d½PI3K > tyrP tyrP > : ¼ k8 ð½IRS2 ½PI3K þ ½IRS1 ½PI3KÞ þ k  8 ð½PI3K_IRSP2  þ ½PI3K_IRSP1 Þ: dt

(4) For the PKR activation, the value of rate constant is borrowed from the activation of PKC in Sedaghat et al. (2002), and the expression for modeling the combinatory effect of JNK and PKR is adapted from Cedersund et al. (2008), which are defined as k12 ¼ k12c þ 10k12c ½JNK P ; k  12 ¼ 0:1  k12c þ 2:99k12c ½AKT p ; where the [JNKP] is normalized to the total amount of JNK. Hence, we have k12

PKR-PKRp ;

k  12

PKR ’ PkRp ;

with 8 d½PKR > P > > < dt ¼  k12 ½PKR þ k  12 ½PKR ; > d½PKRP  > > ¼ k12 ½PKR  k  12 ½PKRP : : dt (5) JNK is activated by MEK, the reaction channels describing the activity are kj5 ;K j5

JNK - JNK P ;

kj6 ;K j6

JNK ’ JNK P ;

with mathematical equations being 8 > kj5 ½JNK½MEK kj6 ½JNK P  d½JNK > > ¼ þ ; > > < dt K j5 þ ½JNK K j6 þ ½JNK P  > kj6 ½JNK P  d½JNK P  kj5 ½JNK½MEK > > ¼  : > > : dt K j5 þ ½JNK K j6 þ ½JNK P  Although there could be different species that are involved in the activation of MAPK/JNK pathways, the detailed kinetic information is unavailable for most of them. Given that we are most interested in the feedbacks and regulations upon IRS1/2, our model simplified the mediators that activate the downstream MAPK/JNK pathway as “component”. Thus, the “Ras” component in the model represents the Ras and Rho GTPase family proteins involving Ras/Rac/Raf, which mediate the MAPK and JNK pathways (Bar-Sagi and Hall, 2000). (6) The dynamics of Fox01 activation through PKR depends on Fox01 itself kf 1

FoxO1-ActiveFoxO1;

kf1

FoxO1 ’ ActiveFoxO1;

with 8 d½FoxO1 > > ¼  kf 1 ½FoxO1 þ k  f 1 ½ActiveFoxO1; < dt d½ActiveFoxO1 > > : ¼ kf 1 ½FoxO1 k  f 1 ½ActiveFoxO1; dt

c

b

[] means concentration of the species.

a

Hargrove (1993), Yao et al. (2008)a k  i1 ¼ 3:0  10  3 210  11 ½ActiveFox01 110  8 þ ½ActiveFox01

ki1 ¼ ki1  k  i1 ½IRS2  2IRS2 IRS2 synthesis (Eq. (33))

Fox01 activation (Eq. (32))

Fox012ActiveFox01

kf 1 ½Fox01  k  f 1 ½ActiveFox01

kf 1 ¼ ð0:1 þ 0:299½PKRP Þk  f 1 k  f 1 ¼ 3

Guo et al. (2006)

C. Huang et al. / Journal of Theoretical Biology 355 (2014) 40–52

where k  f 1 ¼ 3, kf 1 ¼ 0:1k  f 1 þ 0:299k  f 1 ½PKRP .

50

C. Huang et al. / Journal of Theoretical Biology 355 (2014) 40–52

Table 2 Initial concentrations of species. Component Insulin Receptor EN_Receptor IRS1 IRS2 PI3K PI45P2 PI34P2 PI345P3 AKT PKR PKRP Shc GS RasGDP MEK JNK FoxO1

Table 4 Sensitivity analysis of initial conditions.

Initial concentration

∂ ln InRP1 ∂ ln InRP1ð0Þ

References

Based on experimental treatment 9 nM Sedaghat et al. (2002)b 1 nM Sedaghat et al. (2002)b 10 n QUOTE M Sedaghat et al. (2002)b 10/3 nM 1 nM Hatakeyama et al. (2003)b 99.4% Sedaghat et al. (2002) 0.29% Sedaghat et al. (2002) 0.31% Sedaghat et al. (2002) 100% Hatakeyama et al. (2003) 70% 30% 1000 nM Hatakeyama et al. (2003) 10 nM Hatakeyama et al. (2003) 120 nM Hatakeyama et al. (2003)a Raf 120 nM Hatakeyama et al. (2003) 1000 nM Hatakeyama et al. (2003)a ERK 1000 nM

We enhance the initial concentrations of some of the species by Sedaghat et al. (2002) such that it is compatible with that of Hatakeyama et al. (2003). a No evidence or estimations currently available, the values are borrowed from a different system and the original paper is referenced. b Enhanced 10,000–100,000 fold from the original paper referred.

Table 3 Sensitivity analysis of reaction parameters. ∂ ln JNK ¼  1:5684 ∂ ln kj5 ∂ ln RasGDP ¼  1:0019 ∂ ln Vj2 ∂ ln RShGS ¼ 1:0000 ∂ ln k15

∂ ln JNK ¼ 1:5669 ∂ ln Kj5 ∂ ln RasGDP ¼ 1:0000 ∂ ln Kj2

∂ ln RasGDP ∂ ln kj1 ∂ ln GS ∂ ln k17

¼ 1:0019

¼  1:0001

Normalized sensitivity coefficients are used.

(7) FoxO1 activates the transcription of IRS2, and the IRS2 level will be down-regulated if FoxO1 is not activated such that k

i1 φIRS2 ;

IRStyr 1

∂ ln ∂ ln InRP1ð0Þ ∂ ln RasGDP ∂ ln InRP1ð0Þ

where ki1 ¼

2  10  11 ½ActiveFoxO1 1  10  8 þ½ActiveFoxO1

;

k  i1 ¼ 0:003

Since the reaction occurs 15–20 min later, we set a threshold (2.5  10  7 M) on the concentration of active FoxO1 such that it fits the experimental results. The threshold is used to separate the time-scale of transcription (longer time) from the time-scale of other equations modeling signaling processes. All of our numerical simulations are performed using the forward Euler method. A sufficiently small step size of 5  10  4 min is chosen to guarantee convergence of the results. To observe the relations among PKR, IRS1, and IRS2, we plot their time courses in Fig. 5. 3.3. Sensitivity analysis Sensitivity analysis evaluates how the output of a model changes upon small perturbations in the input (parameters or initial conditions). Robust models are not sensitive to small

¼  3:1266

∂ ln PKR ¼ ∂ ln InRP1ð0Þ

¼  3:1003

∂ ln RasGDP ∂ ln IRStyr ð0Þ 1

dXðtÞ ¼ Fðt; XðtÞ; θÞ; dt

3:1420 3:0966

∂ ln IRStyr 1

∂ ln IRStyr ð0Þ 1

¼ 2:3043

∂ ln IRStyr 2 ¼ ∂ ln InRP1ð0Þ

 2:6841

¼ 2:2849

Xð0Þ ¼ X 0 ;

where XðtÞ is a vector of concentrations or percentage of reactants representing the current state of the system, and θ is a vector of parameters. Differentiating both sides of the equation with respect to θ, we have dX θ ∂F ¼ F θ ðt; XðtÞ; θÞ þ ðt; XðtÞ; θÞX θ ðt; XðtÞ; θÞ: ∂X dt Hence the above equation, together with the original system, gives a set of sensitivity equations 8 dXðtÞ > > ¼ Fðt; XðtÞ; θÞ; > > dt > > > < dX θ ∂F ¼ F θ ðt; XðtÞ; θÞ þ ðt; XðtÞ; θÞX θ ðt; XðtÞ; θÞ; ∂X dt > > > > Xð0Þ ¼ X 0 ; > > > : X ð0; θÞ ¼ 0: θ

When the above sensitivity equations are used to compare the relative effect of a parameter change on two or more variables, the sensitivity coefficients must have the same physical dimension or be dimensionless (Turänyi, 1990). Here, we use the normalized sensitivity coefficients. One more step is required to obtain the matrix S of normalized sensitivity coefficients. Sij ¼

ð IRS2 synthesisÞ;

P

changes in inputs, i.e., small changes in the input results in small or no change in the output. Throughout this section, we rewrite the ODE system as follows:

has d½IRS2  ¼ ki1 k  i1 ½IRS2 ; dt

∂ ln IRS1 ¼ ∂ ln InRP1ð0Þ

Normalized sensitivity coefficients are used. Note: InRP1 – phosphorylated once-bound intracellular receptors; IRStyr 1 – tyrosine phosphorylated IRS1 ; PKRP – phosphorylated PKR; IRStyr 2 – tyrosine phosphorylated IRS2 .

k

 i1 φ’ IRS2 ;

¼ 3:1472

θj

X ðt; θÞ ¼ X i ðt; θÞ θ

∂ ln X i ðt; θÞ ∂ ln θj

The most significant sensitivity (relatively large absolute value) coefficients for our model in the time interval [0,10(min)] with constant insulin input are (the absolute value of others are less than 1):Similarly, the sensitivity equations for the initial conditions are 8 dXðtÞ > > ¼ Fðt; XðtÞ; λÞ; > > dt > > > < dX ∂F λ ðt; XðtÞ; λÞX λ ðt; XðtÞ; λÞ; ¼ ∂X dt > > > > Xð0Þ ¼ X 0 ; > > > : X ð0; λÞ ¼ 1: λ

where λ is a vector of initial values. To obtain the normalized sensitivity coefficients for the initial values, we need to perform the normalization as well. The most significant coefficients (with modulus great than 2) are reported in Table 4.

4. Discussion The insulin signaling pathways that control many basic cellular activities are connected to a large number of components. When bound by insulin, the insulin receptor (IR) undergoes receptor auto-phosphorylation to activate its tyrosine kinase activity, which then phosphorylates IRS. At least four of the IR substrates belong to

C. Huang et al. / Journal of Theoretical Biology 355 (2014) 40–52

the IRS group, with IRS1 and IRS2 being predominant and expressed in most tissues (Yang et al., 2010). Tyrosine phosphorylated IRSs are the activators of the phosphoinositide-3 kinase (PI3K), which can catalyze the formation of the lipid second messenger phosphatidylinositol-3,4,5-triphosphoate (PIP3). PIP3 can initiate the downstream AKT and PKC pathway to regulate glucose and lipid metabolism, while PKR is an important downstream effector that can be inhibited by AKT. We showed in previous experimental study that the activated PKR up-regulated the phosphorylations of the serine residues on IRS so as to impair the formation of the activated IRS with tyrosine phosphorylation and thereby inhibit insulin signaling. Furthermore, through binding with Shc, tyrosine phosphorylated IRS also initiates the MAPK pathway, one of whose components is JNK, which in turn downregulates PKR by promoting the inhibitory serine phosphorylation of IRS1 (Yang et al., 2010). Since the MAPK and PKR–IRS pathways are very important they were included in our mathematical model on insulin signaling. Our discrete dynamic model includes these pathways, feedbacks and cross-talks, and shows different physiological functions of IRS1 and IRS2 that emerge from regulatory relationships within the network. The discrete model can be used to determine whether ON/OFF of a component could affect insulin response, based on the logics/wiring in the network. Despite the finding obtained with the discrete model, it is difficult to provide further predictions with the model due to the arbitrarily defined “time steps” and “delays” in the system. The advantage of kinetic model is its ability to provide realistic dynamic profiles and estimate how the output (responses) could be altered by a small perturbation on reaction kinetics. The experimental/clinical intervention of signaling network involves drugs or treatments targeting the activities or bindings of enzymes or kinases, which essentially perturbs the reaction kinetics. Thus we built a kinetic model of insulin signaling based upon the work of Quon (Sedaghat et al., 2002), Hatakeyama (Hatakeyama et al., 2003), and Chan (Yang et al., 2010), that incorporated both the MAPK and PKR–IRS pathways to elucidate not only the regulatory relationship within the network but helped identify points of cross-talks between the signaling pathways that could be potential therapeutic targets. The first part of our model mainly comprises of the first 17 equations in Quon's model, which delineates the binding process of insulin to its receptor, IRS1 phosphorylation and its activation of PI3-kinase, the conversion between PI(3,4,5)P3, PI(4,5)P2 and PI(3,4)P2, and activation of downstream kinase AKT. Hence, our model can reproduce similar profiles as in Quon's study, e.g., Fig. 8 in Sedaghat et al. (2002), but at the same time incorporates IRS2 signaling (see Section 3). Further, our model includes the MAPK model established by Hatakeyama et al. with some parameters adopted from Kholodenko (2000). Note that the unit in Hatakeyama's paper is different from the one in Quon's. Therefore, we apply identities1 nM=s ¼ 60  10  9 M= min and 1=nM=s ¼ 60  109 =M= min to transform Hatakeyama's parameters to fit Quon's model. Moreover, we adjusted some of the parameters in Quon's model to enable our model to capture the dynamics of PKR and IRS that have been shown in previous experimental studies (Wu et al., 2009; Yang et al., 2010). The kinetic model is validated by comparing with the experimental results obtained over a short time-scale (in minutes) of IRS1/2 upon insulin treatment (Wu et al., 2009; Yang et al., 2010) (e.g., compare Fig. 5A and C with Fig. 3 in Wu et al. (2009)). The model is then used to predict the responses at a longer time-scale (several hours) to capture the physiological states (fasting and after feeding). Figs. 2 and 4 show that there is a short-term increase of IRS2 activity right after feeding, but on the long term the activity is down-regulated after feeding and is up-regulated in the fasting state while an opposite trend is shown for IRS1 activity. The kinetic model confirms the different functional roles of IRS1 and IRS2

51

mediated insulin signaling, and provides more details of the dynamics of the system. It could also be used to help identify therapeutic targets for the treatment of insulin resistance through a sensitivity analysis. Despite the aforementioned findings (which are supported by the literature), our kinetic (continuous) model has several disadvantages. First, it requires detailed kinetic parameters and initial concentrations or percentages for each reacting species of the network. Furthermore, some of the parameters are fitted to the experimental data without any means of experimentally validating their values. Second, chemical reactions by nature are stochastic, in particular for those involving reacting species at low concentrations. Our model is in the form of a system of deterministic ODEs, with no noise being introduced. Stochastic modeling captures the detailed reacting process of each molecule. The celebrated Gillespie Algorithm (Gillespie, 2001), which provides an “exact result” in the sense that the probability of a particular realization being computed is the corresponding probability given by the chemical master equation, is every often extremely time-consuming. This is especially true when the whole network is multi-scale, i.e., stiff, which leads to the requirement of tiny step sizes in the algorithm to obtain a solution. Another popular stochastic model is the Langevin equation, which directly introduces Brownian noise into the deterministic ODEs. However, it approximates the chemical master equation well only if concentrations of reacting species are high, which is not the case for our network. From our simulation results, concentrations of reactants vary from 1 nM to 10,000 nM. Compared with the kinetic model, the major advantage of the discrete model is that only qualitative information is required. Thus, it can model and simulate larger networks, given the abundant qualitative interaction maps of biological systems available. Second, the discrete dynamic modeling represents a higher abstraction of biological network, and perturbations on the structure of the network are more easily simulated. The limitations of the discrete model lie in the following aspects (Wu et al., 2009). First, it does not represent or capture the exact timing, but the resulting dynamic patterns provide a dynamic profile evolving with “time steps”, which may contain artifacts due to the network reconstruction and the updating rules at each time step. Second, the discrete models are based solely on the qualitative relationships of the network elements, and thus may lose certain subtle dynamics. It has been shown that in some cases certain kinetic parameters are essential for a system's behavior, which cannot be captured by the qualitative models, e.g., the reaction constants of the phosphatases determining the response time and duration under weak stimulation of the MAPK kinase cascade (Heinrich et al., 2002). For our case, the discrete model provides an overview of the entire system, based upon which, we refine the model with kinetic information, resulting in the detailed kinetic equation system that we can use to study the reaction dynamics and predict potential therapeutic targets. In conclusion, this study constructed systematic models of the insulin signaling network mediated by IRS1 and IRS2, by integrating current knowledge and mathematical models of insulin signaling pathways. Our simulation recapitulated the experimental findings, and suggests that different physiological functions of IRS1 and IRS2 are determined by the wiring and kinetics of the network.

Acknowledgment The research by Huang, Du and Liu is supported by NSF-DMS 084506.

52

C. Huang et al. / Journal of Theoretical Biology 355 (2014) 40–52

Appendix A. Supplementary information Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.jtbi.2014.03.030. References Bar-Sagi, D., Hall, A., 2000. Ras and Rho GTPases: a family reunion. Cell 103, 227–238. Bennett, B.L., Satoh, Y., Lewis, A.J., 2003. JNK: a new therapeutic target for diabetes. Curr. Opin. Pharmacol. 3, 420–425. Brunet, A., Bonni, A., Zigmond, M.J., Lin, M.Z., Juo, P., Hu, L.S., Anderson, M.J., Arden, K.C., Blenis, J., Greenberg, M.E., 1999. Akt promotes cell survival by phosphorylating and inhibiting a Forkhead transcription factor. Cell 96, 857–868. Cedersund, G., Roll, J., Ulfhielm, E., Danielsson, A., Tidefelt, H., 2008. Model-based hypothesis testing of key mechanisms in initial phase of insulin signaling. PLoS Comput. Biol. 4, e1000096. Fritsche, L., Weigert, C., Häring, H.-U., Lehmann, R., 2008. How insulin receptor substrate proteins regulate the metabolic capacity of the liver – implications for health and disease. Curr. Med. Chem. 15, 1316–1329. Gillespie, D., 2001. Approximate accelerated stochastic simulation of chemically reacting systems. J. Chem. Phys. 115, 1716–1732. Gribble, F.M., 2005. Metabolism: a higher power for insulin. Nature 434, 965–966. Gual, P., Le Marchand-Brustel, Y., Tanti, J.-F., 2005. Positive and negative regulation of insulin signaling through IRS-1 phosphorylation. Biochimie 87, 99–109. Guo, S., Dunn, S.L., White, M.F., 2006. The reciprocal stability of Fox01 and IRS2 creates a regulatory circuit that controls insulin signaling. Mol. Endocrinol. 20, 3389–3399. Haeusler, R.A., Accili, D., 2008. The double life of Irs. Cell Metab. 8, 7–9. Hargrove, J.L., 1993. Microcomputer-assistant kinetic modeling of mammalian gene expression. FASEB J.: Off. Publ. Fed. Am. Soc. Exp. Biol. 7, 1163–1170. Hatakeyama, M., Kimura, S., Naka, T., Kawasaki, T., Yumoto, N., Ichikawa, M., Kim, J.-H., Saito, K., Saeki, M., Shirouzu, M., Yokoyama, S., Konagaya, A., 2003. A computational model on the modulation of mitogen-activated protein kinase (MAPK) and Akt pathways in heregulin-induced ErbB signalling. Biochem. J. 373, 451–463. Heinrich, R., Neel, B.G., Rapoport, T.A., 2002. Mathematical models of protein kinase signal transduction. Mol. Cell 9, 957–970. Hirashima, Y., Tsuruzoe, K., Kodama, S., Igata, M., Toyonaga, T., Ueki, K., Kahn, C.R., Araki, E., 2003. Insulin down-regulates insulin receptor substrate-2 expression through the phosphatidylinositol 3-kinase/Akt pathway. J. Endocrinol. 179, 253–266. Kholodenko, B.N., 2000. Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades. Eur. J. Biochem. 267, 1583–1588. Kholodenko, B.N., Demin, O.V., Moehren, G., Hoek, J.B., 1999. Quantification of short term signaling by the epidermal growth factor receptor. J. Biol. Chem. 274, 30169–30181. Kubota, N., Kubota, T., Itoh, S., Kumagai, H., Kozono, H., Takamoto, I., Mineyama, T., Ogata, H., Tokuyama, K., Ohsugi, M., Sasako, T., Moroi, M., Sugi, K., Kakuta, S.,

Iwakura, Y., Noda, T., Ohnishi, S., Nagai, R., Tobe, K., Terauchi, Y., et al., 2008. Dynamic functional relay between insulin receptor substrate 1 and 2 in hepatic insulin signaling during fasting and feeding. Cell Metab. 8, 49–64. Michael, M.D., Kulkarni, R.N., Postic, C., Previs, S.F., Shulman, G.I., Magnuson, M.A., Kahn, C.R., 2000. Loss of insulin signaling in hepatocytes leads to severe insulin resistance and progressive hepatic dysfunction. Mol. Cell 6, 87–97. Paz, K., Hemi, R., LeRoith, D., Karasik, A., Elhanany, E., Kanety, H., Zick, Y., 1997. A molecular basis for insulin resistance. Elevated serine/threonine phosphorylation of IRS-1 and IRS-2 inhibits their binding to the juxtamembrane region of the insulin receptor and impairs their ability to undergo insulin-induced tyrosine phosphorylation. J. Biol. Chem. 272, 29911–29918. Pederson, T.M., Kramer, D.L., Rondinone, C.M., 2001. Serine/threonine phosphorylation of IRS-1 triggers its degradation: possible regulation by tyrosine phosphorylation. Diabetes 50, 24–31. Puigserver, P., Rhee, J., Donovan, J., Walkey, C.J., Yoon, J.C., Oriente, F., Kitamura, Y., Altomonte, J., Dong, H., Accili, D., Spiegelman, B.M., 2003. Insulin-regulated hepatic gluconeogenesis through FOXO1–PGC-1alpha interaction. Nature 423, 550–555. Saltiel, A.R., Kahn, C.R., 2001. Insulin signalling and the regulation of glucose and lipid metabolism. Nature 414, 799–806. Sasaoka, T., Kobayashi, M., 2000. The functional role of Shc in insulin signaling is a substrate of the insulin receptor. Endocr. J. 47, 373–381. Schoeberl, B., Eichler-Jonsson, C., Gilles, E.D., Muller, G., 2002. Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors. Nat. Biotechnol. 20, 370–375. Sedaghat, A.R., Sherman, A., Quon, M.J., 2002. A mathematical model of metabolic insulin signaling pathways. Am. J. Physiol. Endocrinol. Metab. 283, E1084–E1101. Sykiotis, G.P., Papavassiliou, A.G., 2001. Serine phosphorylation of insulin receptor substrate-1: a novel target for the reversal of insulin resistance. Mol. Endocrinol. 15, 1864–1869. Taniguchi, C.M., Ueki, K., Kahn, R., 2005. Complementary roles of IRS-1 and IRS-2 in the hepatic regulation of metabolism. J. Clin. Investig. 115, 718–727. Turänyi, T., 1990. Sensitivity analysis of complex kinetic systems. tools and applications. J. Math. Chem. 5, 203248. White, M.F., 2003. Insulin signaling in health and disease. Science 302, 1710–1711. Wu, M., Yang, X., Chan, C., 2009. A dynamic analysis of IRS-PKR signaling in liver cells: a discrete modeling approach. PLoS One 4, e8040. Yang, X., Nath, A., Opperman, M.J., Chan, C., 2010. The double-stranded RNAdependent protein kinase differentially regulates insulin receptor substrates 1 and 2 in HepG2 cells. Mol. Biol. Cell 21, 3449–3458. Yao, G., Lee, T.J., Mori, S., Nevins, J.R., You, L., 2008. A bistable Rb-E2F switch underlies the restriction point. Nat. Cell Biol. 10, 476–482. Zhang, J., Ou, J., Bashmakov, Y., Horton, J.D., Brown, M.S., Goldstein, J.L., 2001. Insulin inhibits transcription of IRS-2 gene in rat liver through an insulin response element (IRE) that resembles IREs of other insulin-repressed genes. Proc. Natl. Acad. Sci. USA 98, 3756–3761. Zick, Y., 2004. Uncoupling insulin signalling by serine/threonine phosphorylation: a molecular basis for insulin resistance. Biochem. Soc. Trans. 32, 812–816. Zick, Y., 2005. Ser/Thr phosphorylation of IRS proteins: a molecular basis for insulin resistance. Sci. STKE 268, pe4.