A systems biology analysis connects insulin receptor signaling with glucose transporter translocation in rat adipocytes

Type 2 diabetes is characterized by insulin resistance, which arises from malfunctions in the intracellular insulin signaling network. Knowledge of the insulin signaling network is fragmented, and because of the complexity of this network, little consensus has emerged for the structure and importance of the different branches of the network. To help overcome this complexity, systems biology mathematical models have been generated for predicting both the activation of the insulin receptor (IR) and the redistribution of glucose transporter 4 (GLUT4) to the plasma membrane. Although the insulin signal transduction between IR and GLUT4 has been thoroughly studied with modeling and time-resolved data in human cells, comparable analyses in cells from commonly used model organisms such as rats and mice are lacking. Here, we combined existing data and models for rat adipocytes with new data collected for the signaling network between IR and GLUT4 to create a model also for their interconnections. To describe all data (>140 data points), the model needed three distinct pathways from IR to GLUT4: (i) via protein kinase B (PKB) and Akt substrate of 160 kDa (AS160), (ii) via an AS160-independent pathway from PKB, and (iii) via an additional pathway from IR, e.g. affecting the membrane constitution. The developed combined model could describe data not used for training the model and was used to generate predictions of the relative contributions of the pathways from IR to translocation of GLUT4. The combined model provides a systems-level understanding of insulin signaling in rat adipocytes, which, when combined with corresponding models for human adipocytes, may contribute to model-based drug development for diabetes.

Type 2 diabetes is characterized by insulin resistance, which arises from malfunctions in the intracellular insulin signaling network. Knowledge of the insulin signaling network is fragmented, and because of the complexity of this network, little consensus has emerged for the structure and importance of the different branches of the network. To help overcome this complexity, systems biology mathematical models have been generated for predicting both the activation of the insulin receptor (IR) and the redistribution of glucose transporter 4 (GLUT4) to the plasma membrane. Although the insulin signal transduction between IR and GLUT4 has been thoroughly studied with modeling and timeresolved data in human cells, comparable analyses in cells from commonly used model organisms such as rats and mice are lacking. Here, we combined existing data and models for rat adipocytes with new data collected for the signaling network between IR and GLUT4 to create a model also for their interconnections. To describe all data (>140 data points), the model needed three distinct pathways from IR to GLUT4: (i) via protein kinase B (PKB) and Akt substrate of 160 kDa (AS160), (ii) via an AS160-independent pathway from PKB, and (iii) via an additional pathway from IR, e.g. affecting the membrane constitution. The developed combined model could describe data not used for training the model and was used to generate predictions of the relative contributions of the pathways from IR to translocation of GLUT4. The combined model provides a systems-level understanding of insulin signaling in rat adipocytes, which, when combined with corresponding models for human adipocytes, may contribute to model-based drug development for diabetes.
Type 2 diabetes (T2D) 3 is one of our most costly and rapidly spreading diseases. At the heart of this disease lies insulin resistance, an inability of target cells to respond to insulin. This resistance is commonly believed to start in the adipose tissue from which it then spreads to the other target organs, such as liver and skeletal muscle (1)(2)(3). When also these organs become insulin-resistant, the pancreatic beta cells need to significantly increase their insulin production. During this hyperinsulinemic condition, which can last for many years, the glucose homeostasis functions normally, and it is only when the beta cells eventually fail that T2D is manifested. Because of this disease etiology, it is important to understand the insulin signaling network in adipocytes because such an understanding may provide ideas for a possible cure for T2D before it has manifested.
The intracellular insulin signaling network is a complex web of thousands of proteins, and even though many questions remain unanswered, several core pathways have been established (for a review, see Ref. 4). All of these pathways originate with binding of insulin to its receptor (the insulin receptor (IR)), and the perhaps most well studied pathway leads to increased glucose uptake via translocation of glucose transporter 4 (GLUT4) to the plasma membrane. In this pathway, binding of insulin to IR leads to autophosphorylation at IR tyrosine residues, which in turn induces tyrosine phosphorylation of insulin receptor substrates (IRSs). This is followed by a sequence of steps, including increased production of membrane-bound phosphatidylinositol 3,4,5-trisphosphate (PIP 3 ), that lead to a transient membrane translocation and phosphorylation of protein kinase B at threonine 308 (PKB-Thr-308) (5). PKB is also phosphorylated by insulin at another site, serine 473 (PKB-Ser-473). Unlike for Thr-308, PKB-Ser-473 phosphorylation is not mediated via IRS but most likely via a pathway involving the mammalian target of rapamycin complex 2 (mTORC2) (6). There may also be an increased autophosphorylation of PKB-Ser-473 caused by PKB-Thr-308. Once activated, PKB further phosphorylates several sites at the downstream target Akt sub-strate of 160 kDa (AS160), which is one of the key factors required for GLUT4 translocation to the plasma membrane (7). In addition to this well established pathway, there are several other factors regulating GLUT4 translocation that are less defined, including not fully elucidated modifications at the membrane (8,9) and priming of the insulin-containing vesicles (Ref. 10; for a recent review, see Ref. 11). Apart from this lack of mechanistic details, the quantitative and relative contribution of all the upstream signaling factors is still not established. In other words, even though many of the intermediates and processes that are involved are more or less identified, how these processes interact, which processes are rate-limiting, and how these processes can be modified in healthy and diseased conditions are still far from fully understood. In short, a systems-level understanding is missing.
A systems-level understanding in the field of biology, also called a systems biology understanding, rests upon two pillars: (i) the simultaneous collections of systems-wide and quantitative experimental data from comparable conditions and (ii) analysis of these data using mechanistic mathematical modeling. Importantly, fully developed mechanistic mathematical modeling goes beyond mere correlations and statistics, which only describe interrelations between data but not the interactions of the underlying biological system. In other words, a mechanistic systems biology model tests and develops mechanistic descriptions of a system. Once an acceptable description has been generated, this can be simulated in a computer to produce in silico experiments, which can be used to e.g. plan, design, and predict future experiments and serve as estimates of non-measurable properties (12,13).
This systems biology approach has been used for several decades to unravel insulin binding to the IR (14,15), the insulin signaling pathway downstream of the IR (16), and GLUT4 dynamics in both primary rat adipocytes (17)(18)(19)(20)(21) and cultured 3T3-L1 cells (22). In human adipocytes, these models have lately been expanded to include IR activation and the entire pathway leading to glucose uptake (16,(23)(24)(25). These recently developed human models are now at use in pharmaceutical drug development, describing both normal signaling based on data collected in adipocytes from healthy subjects and insulinresistant signaling with data collected in adipocytes from diabetic subjects. Also, these human models translate the results based on cellular data to the whole-body level (16,23). In rat adipocytes, in contrast, such interconnective models based on time-resolved data do not exist even though some detailed models exist for certain subsystems. One such model describes the GLUT4 translocation process (Fig. 1, green box). This model has provided non-trivial conclusions regarding the necessary number of pools involved in this GLUT4 recycling machinery based on live-cell image analysis with temporal and spatial resolution, which allows for direct measurements of many of the parameters included in the model (26). For the early insulin signaling events in rat adipocytes, the counterintuitive dynamics seen between IR and IRS1 have also been unraveled using modeling (27) (Fig. 1, blue box), also here providing non-trivial insights regarding the role and quantitative relations between different pools of IRs. However, unlike for human adipocytes, no model for primary rat adipocytes that connects the upper IR-IRS1 dynamics with the ultimate GLUT4 translocation dynamics exists (Fig. 1, middle questions). This is a critical lack because much drug development still is done using rodents, and many key questions only can be tackled by models that describe the entire insulin signaling network, such as (a) which upstream pathways are regulating which steps of the GLUT4 vesicle translocation, (b) what is the quantitative contribution of these upstream pathways to the GLUT4 vesicle translocation, and (c) is our current knowledge sufficient to describe all available data for insulin signaling in rat adipocytes, or are additional not yet elucidated mechanisms necessary?
Here, we combine existing data for IR and IRS1 activation and GLUT4 translocation dynamics with new measurements of the key intermediate steps PKB-Thr-308, PKB-Ser-473, and AS160-Thr-642 phosphorylations in primary rat adipocytes. These combined data then allow us to construct a model, based on comparable rat adipocyte data, for the entire pathway leading from IR activation to GLUT4 translocation. Different versions of acceptable parameter values and model structures are used to quantify the importance of the upstream signals for the different steps in the GLUT4 translocation machinery. The resulting model is a quantum step forward toward a systems-level understanding of insulin signaling also in rat adipocytes, bringing us one step closer to knowledge-driven drug development that bridges between the preclinical and clinical phases.      Figure 1. Overview of the project. Models for the insulin binding and the subsequent activation of IR and IRS1 (blue box) and the dynamic GLUT4 vesicle translocation (green box) already exist. However, the connecting pathways between them are not entirely known, and several questions exist. Which are the pathways? How much does each pathway contribute? Is our mechanistic understanding of the insulin signaling chain sufficient to describe all our observations? We herein develop and test these three submodels (Figs. 2-4), integrate them into a combined model (Fig. 5), test the model with respect to new validation data (Fig. 6), and identify unique predictions to quantify the properties of the three pathways (Fig. 7). The star in the upper model graph represents an unknown type of activation. IRi, internalized insulin receptor; X, unknown protein; YP, known activation.

Module 1, describing the IR-IRS1 dynamics
The new model is now described and compared with data for each of the three Modules 1-3, whereafter the combined model is analyzed. The dynamics of Module 1, activation, internalization, and recycling of the IR and the subsequent activation of IRS1, have been unraveled in both human adipocytes (15) and rat adipocytes (27). Because we herein study insulin signaling in rat adipocytes, the insights and corresponding model from rat adipocytes were used. In these cells, the IR and IRS1 appear to behave in a counterintuitive manner where the peak of phosphorylation of IR (Fig. 2B) occurs later than the peak of phosphorylation of the substrate (IRS1) (Fig. 2C). In Nyman et al. (27), this behavior could be explained with a negative feedback or with the fact that there are different pools of IRs that have different dynamic behaviors. These suggestions were tested in detail in Nyman et al. (27), which also included data on the fraction of IR remaining in the plasma membrane after insulin stimulation (Fig. 2D), and here we chose to include the simplest acceptable explanation: a negative feedback via an unknown protein X ( Fig. 2A). This negative feedback from downstream IRS1 back to IRS1 could be interpreted as (i) a decrease in the kinase activity of a protein that normally increases the activity of IRS1, (ii) an increase in the activity of a phosphatase that normally dephosphorylates IRS1, or (iii) an increase in specific serine phosphorylations (i.e. Ser-636) that are believed to decrease the activity of IRS1 (28). In the first module, we use the model equations and parameters for IR and IRS1 from Nyman et al. (27), and we verified that the model simulations (Fig. 2, B-D, lines) are in good agreement with the data (error bars) ( 2 ϭ 22.02 Ͻ cutoff ϭ 35.17; p Ͻ 0.05) (see Materials and methods in Ref. 12). The cutoff is the inverse of the 2 cumulative distribution function and was computed with the degrees of freedom specified by the number of data points (26) and the chosen probability (95%). The experimental data in Module 1 for IR and IRS1 dynamics were obtained with 100 nM insulin, and we therefore used 100 nM insulin in the model simulations. The data in the different modules are taken from independently conducted experiments with different insulin concentrations, and we have altered the insulin concentration in the model accordingly when simulating the different experiments.

Module 2, signaling via PKB and AS160
In Module 2 ( Fig. 3A), we included known pathways that bridge between Module 1, for the IR-IRS1 module, and Module 3, for GLUT4 translocation. The most well established pathways go via PKB, which is phosphorylated both at site Thr-308, via activation of PIP 3 and phosphoinositide-dependent protein kinase-1, and site Ser-473, via mTORC2. Downstream of PKB, we also included AS160, which is known to be essential for regulation of GLUT4 translocation (7). To be able to find values for the parameters in Module 2, we performed experimental measurements of the dynamic response to insulin of PKB phosphorylation at sites Thr-308 and Ser-473 and AS160 at site Thr-642. Isolated, primary rat adipocytes were stimulated with insulin (28 nM) for 0, 0.5, 1, 2, 3, 5, 10, or 20 min. Cell lysates were subjected to Western blotting, and protein phosphorylation was quantified using phosphospecific antibodies; all samples were normalized to ␤-actin (n ϭ 4 independent experiments) ( Fig. 3, B-D). For this experiment, the insulin concentration was 28 nM, which was also used in the model when fitting the model parameters. A good agreement was found between model simulations and the experimental data of PKB and AS160 phosphorylation (Fig. 3, B-D, lines and dots) ( 2 ϭ 27.98 Ͻ cutoff ϭ 36.42; p Ͻ 0.05) (see Materials and methods in Ref. 12). The 2 cutoff was once again computed with the degrees of freedom specified by the number of data points (24 points) and the chosen probability (95%). The remaining parts of Module 2 are easiest to describe in connection to a description of the inputs to Module 3.

Module 3, connecting with GLUT4 vesicle translocation
In Module 3 (Fig. 4A), we used the dynamic GLUT4 vesicle translocation model developed in Stenkula et al. (26) and tested three input connections from Module 2: 1) from AS160-Rab-GTP, 2) from an AS160-independent pathway downstream of PKB, which we refer to as "PKBd" and which is given by a delayed state downstream of the sum of the two states that are phosphorylated at Thr-308, and 3) from a plasma membranelocated pathway in close proximity to the IR, which we refer to as "PMR." The first connection, from AS160-RabGTP, is described in a simplified view compared with the actual mechanisms, which are only partially understood. The more detailed underlying mechanisms include the following. (i) AS160 in the non-insulin stimulated state, through its GAP domain, maintains the inactive GDP-bound form of Rab. (ii) AS160 upon insulin stimulation becomes phosphorylated (Thr-642), which inhibits its GAP activity. (iii) This inhibited activity allows the active form of Rab (RabGTP) to increase, which leads to association of the GLUT4-containing vesicles with the plasma membrane (7). The net effect of these three facts (i-iii) is that insulin stimulates GLUT4 translocation to the plasma membrane via the intermediates AS160 and Rab. We modeled these steps via a simple cascade describing two inhibitory steps, which thus cancel each other out, giving a net positive effect of insulin on GLUT4 translocation (Figs. 3A and 4A and combined in Fig. 5). More specifically, we allowed the effect of AS160-RabGTP to increase both the rates V1 and V2, which represent the translocation of GLUT4 to the plasma membrane followed by release of GLUT4 monomers (V1) and translocation of GLUT4 storage vesicles (GSVs) leading to formation of GLUT4 clusters in the plasma membrane (V2), as defined in Stenkula et al. (26) (Fig. 4A).
The second connection, from PKBd, describes a PKB-Thr-308 dependent, but AS160-independent, pathway. This connection was included to describe differences between the insulin-induced effects measured for V1 and V2: V1 had a much higher -fold increase with insulin compared with V2 (Fig. 4B). The -fold increase was measured in Stenkula et al. (26) where it was concluded that translocation of GSVs to the plasma membrane followed by release of GLUT4 monomers (V1) is the fusion reaction that responds the most to insulin, whereas translocation of GSVs forming GLUT4 clusters (V2) at the plasma membrane is the most prominent fusion reaction in the non-stimulated state. The measured -fold from basal to peak for V1 is the factor 60, and for V2 the -fold is the factor 2. The simulated -fold values are 58.1 and 1.98, respectively, which thus fit the experimental data well. The sum of V1 and V2 was also measured (Fig. 4C).
The third connection, PMR, was included to illustrate an insulin-induced effect associated with the plasma membrane downstream of the IR but also affecting the measured insulininduced effect at the rate of GLUT4 molecules moving out of GLUT4 clusters into the pool of GLUT4 monomers (Vr) and the measured inhibitory effect at the rate of GLUT4 endocytosed from GLUT4 clusters into endosomes (Ve). These rates were measured in experiments with and without inhibition of endocytosis for different levels of relative exposure in Lizunov et al. (29), as summarized in Fig. 4, E and F, under insulin conditions and with the individual relative exposures in Fig. 4, G-J, under basal conditions. All data in Fig. 4 were used to parametrize the three connections between upstream insulin signaling and GLUT4 translocation. As can be seen in Fig. 4, B-J, the model can describe the data well according to a visual inspection; a 2 test is not possible because we lack uncertainties for some of the data points. A total of 98 data points were used in Module 3 (92 visualized in graphs plus four rate constants and 2-fold constants). In both studies describing experimental data for Module 3 (Fig. 4, B-J, error bars), the cells were incubated with an insulin concentration of 70 nM, which also was used in the simulations.

Verification of the novel combined insulin signaling model
We combined the three modules into a model for insulin signaling from IR activation to complete GLUT4 translocation in rat adipocytes (Fig. 5). The model is in simultaneous agreement with all experimental data used for model development Isolated primary rat adipocytes were incubated with insulin (28 nM) for the indicated time (0 -20 min). Cell lysates were subjected to Western blotting, and protein phosphorylation (p) was detected by phosphospecific antibodies toward AS160-Thr-642, PKB-Ser-473, and PKB-Thr-308. ␤-Actin was used as loading control. n ϭ 4 independent experiments (E). PIP2, phosphatidylinositol 4,5-bisphosphate; PDK1, phosphoinositide-dependent protein kinase-1.  . This is, to our knowledge, the first model of insulin signaling in rat adipocytes that describes the entire dynamic and quantitative pathway from IR activation to GLUT4 translocation.

Systems biology analysis of insulin signaling
To test the soundness of the model, we studied its behavior in new scenarios for which it had not been tested. Because the key new aspect of our model is the upstream signaling feeding into the GLUT4 cycling, we identified data from the same experimental system that perturbed these pathways differently.
These data are based on a study where rats have been treated with dexamethasone prior to isolation of adipocytes and insulin stimulation (30). Dexamethasone treatment inhibits insulinstimulated PKB phosphorylation at Thr-308 but only marginally (if at all) at Ser-473. Dexamethasone also decreases the total expression of PKB (30). In our model, we implemented these changes by reducing the two parameters k308_1 and k308_2 (Fig. 6A) so that PKB-Thr-308 phosphorylation is reduced to 50 or 70% (Fig. 6B, dashed line), which is the experimentally observed inhibition range (Fig. 6B, error bar). We also reduced the total amount of PKB in the model according to the range in data (37-51% reduction). We used these changes in PKB total amount and Thr-308 phosphorylation to predict the resulting alterations in glucose uptake and plasma membrane localization of GLUT4 for the same parameter values that were used for the plots in Figs. 2-4. Glucose uptake is assumed to be proportional to the total amount of GLUT4 in the plasma membrane. For these chosen parameter values, the GLUT4 plasma membrane pool is reduced by 16 -22% with dexamethasone inhibition (Fig. 6C). To account for the uncertainties in all parameters and data points, we used the method proposed in Cedersund (13). The boundaries for the prediction are identified with a direct optimization where the reduction in GLUT4 in the plasma membrane/glucose uptake is stepwise changed while all parameters are optimized for the lowest possible cost, i.e. to find the best agreement with all data (Fig.  6D, line) ("Materials and methods"). We chose 105 as the cutoff level for the cost, and this gave an interval of 0.15 and 0.26, i.e. a 15-26% reduction (Fig. 6D, blue area). The predicted inhibition of GLUT4 in the plasma membrane is thus well determined and can be compared with the corresponding experimental value (30) (Fig. 6D, red area). The predicted value overlaps with the experimental value (Fig. 6E), which means that the model is in good agreement with data not used for fitting the parameters. The insulin concentration used when acquiring the experimental data was 7 nM, or 1000 microunits/ml, and thus the same concentration was used to simulate these model predictions.

Predicting the upstream effectors to GLUT4 translocation
Having established some faith in the model, we therefore moved on to characterize some predictions that are of high biological interest for the new aspects of the model: the three input signals that feed into the GLUT4 translocation machinery. Using our model, we searched for new parameters that would "stretch" the prediction of the -fold increase of the three input signals, RabGTP, PKBd, and PMR, in response to 70 nM insulin. For that, we used the same method as above (13) to stepwise perform the stretch and at the same time stay in good agreement with the experimental data. This resulted in a 15-28-fold increase from basal to insulin steady state of PKBd (Fig. 7A), an interval of 1.4 -2.25 for PMR (Fig. 7B), and an interval of 2.65-3.35 for RabGTP (Fig. 7C).
In summary, our analysis shows that the contributions of the upstream signals to GLUT4 vesicle translocation, RabGTP, PKBd, and PMR, all are constrained by available data. None of these three upstream signals have been measured in experiments. Therefore, these new predictions provide a first assessment of the dynamic range of these three signals in response to insulin stimulation.

Discussion
Herein, we have developed the first mathematical model based on time-resolved data that describes the entire insulin signaling pathway from IR activation to GLUT4 exocytosis in rat adipocytes. The model was developed in three modules: Module 1, IR-IRS dynamics (Fig. 2); Module 2, PKB dynamics and the pathways that are involved in GLUT4 translocation (Fig. 3); and Module 3, containing the actual GLUT4 translocation dynamics (Fig. 4). The combined model (Fig. 5) is developed based on Ͼ140 data points, including immunoblotting on the IR level and live-cell imaging data of GLUT4 exocytosis from previous studies, as well as new data generated to obtain temporal resolutions of several insulin signaling intermediates (Fig. 3, B-D). Our model can describe not only all of these data simultaneously (Figs. 2, B-D; 3, B-D; and 4, B-J) but also independent validation data without parameter refitting (Fig. 6, D . C, this PKB-Thr-308 reduction together with the reduction of total PKB in turn leads to a predicted decline in GLUT4 vesicle translocation to the plasma membranes by 16 -21% at insulin steady state compared with adipose cells treated with insulin alone for one set of parameters. The range for this one parameter was obtained by simulating either with the maximal S.E. boundaries (64% reduced PKB expression, 50% reduced PKB-Thr-308 activation; yellow dashed line) or minimal S.E. boundaries (49% reduced PKB expression, 30% reduced PKB-Thr-308 activation; red dashed line). D, to find the overall maximal and minimal values of the predicted GLUT4 reduction (blue interval), the prediction was lowered in steps while optimizing all parameters to data to see how long this model structure can agree with estimation data (i.e. how long the blue curve can stay below the cutoff) (13). Assuming that the fraction of GLUT4 in the plasma membrane is proportional to the glucose uptake, this prediction (15-26%) overlaps almost completely with the experimentally observed reduction (red interval; 17-50%). E, alternative visualization of the experimental validation in D.  (Fig. 7). In summary, the combined model provides a new level of understanding for insulin signaling in rat adipocytes: a quantitative systems-level understanding that both summarizes available knowledge and data and is able to predict new experiments that come as consequences of the new systems-level understanding. The work with defining the connecting pathways from the insulin signaling IR-IRS1 module down to the GLUT4 translocation module has lead us down to three independent input signals, and all of these input signals are based on previously proposed mechanistic hypotheses and literature. In other words, none of the three linking mechanisms are new hypotheses. One input signal is centered around AS160. The identification and characterization of AS160 as a key regulatory mediator downstream of PKB was a major discovery in understanding the mechanisms of insulinstimulated glucose transport in adipocytes (7). Knockdown of AS160 increases the pool of ready-to-fuse vesicles, hence shifting the GLUT4 equilibrium to the plasma membrane already in the non-stimulated state. However, even with diminished AS160 levels, insulin has been shown to increase GLUT4 exocytosis (31), and reports have suggested that another PKB substrate, independent of AS160, exists (32). This is the basis for the presence of the second input signal, an AS160-independent pathway PKBd. In the model, PKBd directly feeds into the GSV fast fusion events tightly regulated by insulin, which fits the idea of PKB as a prominent regulator of GSV fusion close to the plasma membrane. The fact that our model supported PKB-Thr-308 activation as essential for downstream signaling leading to GLUT4 translocation is in line with previous reports from rat adipocytes (33) and herein served as a verification of the model design. The third and last input signal is centered around plasma membrane activation: PMR. This is based, in part, on an elegant study by Koumanov et al. (8) where the rate-limiting step of insulin-regulated GLUT4 exocytosis involved a plasma membrane activation, which was supported in later studies (26). This third signal could potentially involve clustering of SNARE proteins at the plasma membrane (for a review, see Ref. 34), which could increase the probability of fusion, similar to exocytosis of insulin granules (35). Also, previous findings in cultured adipocytes suggest a direct interaction of the IR and the SNARE protein Munc18c required for complete GLUT4 exocytosis (36) as well as local PIP 3 accumulation following insulin stimulation (37). In our model, PMR was included as a plasma membrane-associated step downstream of the IR that represents the insulin-induced effects associated with the plasma membrane, which in the model inhibits the rate of GLUT4 endocytosis, increases the rate of monomer formation from clusters, and increases the rate of GLUT4 exocytosis into monomers (Fig. 5).
Although the presence of these three input signals is based on existing literature, we provide a first systems-level description of how they could act together. In this model, the three input signals act in different ways to stimulate different parts of the GLUT4 subsystem. In most cases, a single input signal is responsible for the complete change observed in a specific rate: PMR is alone responsible for the change observed in endocytosis (Ve) and cluster-to-monomer breakdown in the membrane (Vr), and RabGTP is alone responsible for the increased rate of GLUT4 cluster formation (V2) (Fig. 4A). In contrast, for the rate of monomer formation from exocytosis (V1), all three To predict the -fold values of the three inputs to the GLUT4 translocation module, a stepwise increase/reduction of the -fold values was performed in the same manner as in Fig. 6D (13). The -fold values are forced to attain certain values using an extended cost function as explained under "Material and methods." As long as the -fold is acceptable (final cost lower than the predetermined cutoff), the -fold is further increased or reduced, which eventually creates a range of acceptable -fold values. Finite ranges for the inputs, PKBd (A), PMR (B), and RabGTP (C), could be obtained, implying that all three predictions are identifiable for the given data and model structure.
input signals contribute. Nevertheless, their contributions are different because their dynamic range is different. The experimentally observed dynamic range (-fold change from basal to insulin steady state) for V1 is high: ϳ20 (26). This dynamic range is higher than what PMR and RabGTP can have because they are constrained by the data in Fig. 4, E-J, (29) and by the data for V2 (Fig. 4B, purple), respectively. In contrast, the dynamic range for PKBd goes much higher because PKB-Thr-308 has a large dynamic range (Fig. 3B) and because PKBd needs to account for the difference between the large dynamic range seen for V1 and the small dynamic ranges required for RabGTP and PMR. In other words, neither RabGTP nor PMR alone or their combination could act to produce the large dynamic range seen in V1, but PKBd must be present as well. This reasoning based on the data alone illustrates how the predictions produced in Fig. 7 can be understood intuitively even though a model-based analysis is needed to ensure that the reasoning holds true also when taking all quantitative features in all available data into account. In summary, we can thus answer the three questions in Fig. 1: our mechanistic understanding is sufficient, there are three connecting pathways (Fig.  5), and their contribution is specified by how they enter the GLUT4 module (Fig. 5) and by their -fold responses (Fig. 7).
Here, we have focused on data for primary rat adipocytes. It is likely that the distinct architecture of these primary cells, with a thin cytosolic rim surrounding the lipid droplet, in contrast to cultured adipocyte cell lines where cells contain multiple lipid droplets and an increased cytosolic domain implies more or less different insulin signaling networks leading to GSV exocytosis. For example, it is well known that 3T3-L1 cells have a substantial amount of glucose transporter 1 (GLUT1) (38), which is not regulated by insulin but does contribute to glucose uptake. In this study, we have therefore concentrated our analysis exclusively on experimental data from primary rat adipocytes and not built further on models developed for 3T3-L1 cells (22,39). In general, we believe that modeling of this kind requires that all data have been collected from the same cell type under comparable conditions where all existing differences between the experiments can be described within the model. For the same reason, we do not herein make use of any of the similar data and model results available for human adipocytes (16).
As with all models, the model developed herein is not a final or complete description of the studied system but is a simplified description based on a limited set of experiments. In other words, this model only includes such mechanisms that are found necessary to describe the current set of data and not all known mechanisms and variations. This choice of a minimal acceptable model should minimize the prediction uncertainty, at least for those predictions that are of the same character as the already existing data (12). Our model-based analysis moves the understanding of insulin signaling in primary rat adipocytes to a new level: to a quantitative systems-level understanding. Our analysis is systems-level because we (i) collect data from all levels of insulin signaling (early receptor activation, downstream signaling, and GLUT4 translocation) and (ii) describe all of these data simultaneously using a single mathematical model. Our analysis is quantitative because the model describes quantitative features in the data, such as correct -fold changes and the detailed shapes of the responses. To the best of our knowledge, this is the first instance of such a model for this system, and this new level of understanding opens new possibilities. We can therefore now do useful predictions with corresponding uncertainties, which is the model's propagation of all of the uncertainties in the original experimental data to this particular prediction.
The model developed herein for insulin signaling in rat adipocytes can also be used together with models developed for the same kind of data from human adipocytes (such as Refs. 16 and 25). These human models exist in an insulin-resistant/diabetic version and can thus be directly used to study which perturbations drive the model from a diabetic state to a normal insulin response. For the model developed herein, we need corresponding data from adipocytes from rats that display a diabetic phenotype, for example Zucker diabetic fatty rats, that are commonly used in drug development. With this additional information, we can use the model to predict drug responses in both human adipocytes from diabetic patients and in rat models of diabetes and search for drugs that have an effect in both these systems. This is important because the drug development pipeline includes preclinical tests in rodents where an effect must be shown before the drug is taken to clinical trials.

Mathematical modeling methods
We have used ordinary differential equations when formulating the model, and the forms of the equations are directly given by the interaction graphs in Figs. 2A, 3A, 4A, and 5. This means that the IR-IRS module has the following equations.
where IR, IRp, IRi, IRS, IRSp, X, and Xp are states with initial conditions as given, and v1a, v1b, vR, v2, v3, vm3, v4, and vm4 are reaction rates given by the following equations. where ins is the input (insulin added to cells), and ik1, ik1basal, ik2, ikR, ik3, ikf, ikm3, ik4, and ikm4 are model parameters. The output of the model is compared with data. For the IR-IRS module, the following output was measured.
where measIRS, measIR, and measIRmem are variables that are compared with the data in Fig. 2B, and scaleIRS, scaleIR, and scaleIRmem are parameters used to scale between model simulations and data.
The values of the rate constants are hard to measure experimentally. Instead, these values are estimated by finding the values of the parameters that maximize the agreement between simulated outputs of the ordinary differential equation model and the experimental data. This is done by minimizing the cost function, V(p), that sums the squares of the residuals between data and model simulations weighted by the data variance.
where t denotes measured time points, and p denotes model parameters. Additionally, y is the experimental data, ŷ is the simulated variable value, and S.E. is the standard error of the mean of data. The sum i is over all measured variables for all measured time points. To statistically test the quality of the model's agreement with the data for a parameter p, we use a 2 test: the cost function V(p) gives us the 2 value, which is compared with an inverse of the 2 cumulative distribution function computed using the degrees of freedom given by the number of data points (12). In this report, Module 1 and Module 2 were statistically evaluated. Module 3 was not eligible for statistical evaluation because several data points are missing S.E. values. The following is one set of acceptable values for the parameters of the IR-IRS module.
To simulate the model for a given experiment, the input (ins) is changed to the value that corresponds to the added concentration of insulin in that particular experiment. Insulin concentrations in the data used for estimation are in the range 7-100 nM. For the prediction of the -fold increase of the three branches of insulin signaling (Fig. 7), we chose to use 70 nM insulin, which is the insulin concentration used in the experiments of Module 3 (70 nM). In such predictions, we account for the fact that we do not have unique values of the parameters and that some parameters even may be unidentifiable from the available data. This is done by using the method in Cedersund (13) and especially the version advocated in Cedersund et al. (40). In this version, one adds an extra term to the cost function to create an extended cost function.
where w is a large number, h(p) is the studied model property (e.g. the -fold increase), and h is the value that we want the model property to obtain. If w is large enough, a minimization of Ṽ ensures that h(p) ϭ h and that one seeks the best agreement to the estimation data for this given value of the model property. The agreement with the data, V(p), is plotted as the y axis in Figs. 6D and 7, A-C; h(p) is the x axis. All mathematical analyzes, i.e. model simulations and parameter estimations, have been carried out in Matlab together with the SBtoolbox2 package (41, 42). The scripts for these analyses and how to simulate all the figures herein can be found in the supplemental material.
Cell preparation and Western blotting-Adipose cells were isolated from epididymal tissue as described previously (43). Following isolation, cells were suspended (10% suspension) in Krebs-Ringer buffer containing 25 mM Hepes, pH 7.4, 200 nM adenosine, and 1% BSA (w/v) and stimulated with insulin as indicated in the figures (at 37°C with shaking at 150 cycles/ min). To stop incubations, cells were washed in Krebs-Ringer without BSA on ice and subsequently lysed in a buffer containing 50 mM Tris/HCl, pH 7.5, 1 mM EGTA, 1 mM EDTA, 1 mM sodium orthovanadate, 10 mM sodium ␤-glycerophosphate, 50 mM sodium fluoride, 5 mM sodium pyrophosphate, 0.27 M sucrose, 1% Nonidet P-40, 1 mM dithiothreitol (DTT), and Complete protease inhibitor mixture (one tablet/50 ml). Lysates were centrifuged for 15 min at 1000 ϫ g, and protein concentrations were determined by the method of Bradford (44). Protein (10 g/sample) was heated at 95°C for 2 min in SDS sample buffer and subjected to polyacrylamide gel electrophoresis on precast Bio-Rad gradient gels and electrotransfer to nitrocellulose membrane. Membranes were blocked for 30 min in 50 mM Tris/HCl, pH 7.6, 137 mM NaCl, and 0.1% (w/v) Tween 20 (TBS-T) containing 10% (w/v) milk powder. The membranes were then probed with the indicated antibodies in TBS-T containing 5% (w/v) milk powder or 5% (w/v) BSA for 16 h at 4°C. Detection was performed using horseradish peroxidase-conjugated secondary antibodies and the enhanced chemiluminescence reagent. The signal was visualized using a Bio-Rad imaging camera, and band intensities quantified using Bio-Rad imaging software.

Systems biology analysis of insulin signaling
Author contributions-G. C., K. G. S., N. B., and E. N. contributed to study design, researched data, discussed and interpreted data, and wrote the manuscript. All authors reviewed the results and approved the final version of the manuscript.