Continuous variable responses and signal gating form kinetic bases for pulsatile insulin signaling and emergence of resistance.

Understanding kinetic control of biological processes is as important as identifying components that constitute pathways. Insulin signaling is central for almost all metazoans, and its perturbations are associated with various developmental disorders, metabolic diseases, and aging. While temporal phosphorylation changes and kinetic constants have provided some insights, constant or variable parameters that establish and maintain signal topology are poorly understood. Here, we report kinetic parameters that encode insulin concentration and nutrient-dependent flow of information using iterative experimental and mathematical simulation-based approaches. Our results illustrate how dynamics of distinct phosphorylation events collectively contribute to selective kinetic gating of signals and maximum connectivity of the signaling cascade under normo-insulinemic but not hyper-insulinemic states. In addition to identifying parameters that provide predictive value for maintaining the balance between metabolic and growth-factor arms, we posit a kinetic basis for the emergence of insulin resistance. Given that pulsatile insulin secretion during a fasted state precedes a fed response, our findings reveal rewiring of insulin signaling akin to memory and anticipation, which was hitherto unknown. Striking disparate temporal behavior of key phosphorylation events that destroy the topology under hyper-insulinemic states underscores the importance of unraveling regulatory components that act as bandwidth filters. In conclusion, besides providing fundamental insights, our study will help in identifying therapeutic strategies that conserve coupling between metabolic and growth-factor arms, which is lost in diseases and conditions of hyper-insulinemia.

Signaling cascades are essential for regulating cellular processes, and decades of work have unraveled molecular and biochemical mechanisms that constitute them. However, kinetic parameters that define emergent properties of signaling networks and therefore predict regulatory nodes are poorly understood. While independent experimental and mathematical approaches have provided valuable insights (1–6), studies that capture dynamics and complexities of signaling architecture vis-à-vis physiological variations in input strengths are far fewer. Not only would these reveal fundamental kinetic considerations that determine signal topology but also inform about reactions/components that could emerge as therapeutic targets.

Evolutionarily conserved insulin signaling (IS) is essential for cellular/organismal metabolism and growth (7–9). Aberrant IS is associated both causally and consequentially with growth abnormalities, inflammation, accelerated aging, and diseases including metabolic disorders and cancer (10–13). Genetic perturbations and omics-based studies have elucidated importance of key phosphorylation events in response to insulin stimulation (14–16). Recent reports have provided crucial insights into physical protein interactome, temporal changes in phospho-proteome, and kinetic constants, viz. half-maximal time to reach peak response and half-maximal effective concentration (EC50) (17, 18). However, kinetic parameters that govern network properties of IS as a function of normo-insulinemic and hyper-insulinemic states, which could collectively determine physiological and pathophysiological outcomes, is still lacking.

Our current understanding largely stems from studies, which have used either supraphysiological or static concentrations of insulin (19–21). This is in contrast to the physiological setting wherein circulating insulin concentrations vary drastically from being low-pulsatile (∼0.1 nM) to high-biphasic (∼1.0 nM) in fasted and fed states, respectively (22). Moreover, kinetic criteria that either encode fasted-to-fed transitions or drive pathological manifestations of IS, as in diabetes and obesity (23, 24), are unknown. Since IS can be broadly divided into metabolic and growth-factor arms (25), if/how the flow of information is stratified and maintained under normal and hyper-insulinemic states remains to be unraveled.

Mathematical approaches to model cellular signaling have gained traction in the recent past to understand the dynamics and also to provide predictive parameters that define topology of signaling network (26–29). Earlier such attempts to determine kinetics of IS have largely employed “averaged” measures to define the behavior of the system (3, 17). Notably, given the fluctuations in insulin levels and inherent noise in signaling, there are no reports that have computed kinetic parameters, which capture emergent properties of IS. Specifically, while there have been simulation-based approaches to define dose-to-duration effects and kinetic insulation on synthetic signaling networks (30), such principles have not been applied to complex biological cascades.

In this current study, we unravel network properties of IS under physiological concentrations of insulin and reveal kinetic basis for emergence of memory and resistance. Our study utilizes parameters such as kinetic barriers and connectedness in the network to address how signaling topology is maintained. Notably, we describe the importance of dynamic range and pulsatility in signaling, which generates memory as well as couples the metabolic and growth-factor arms.

Results

Distinct Kinetics of Signaling in Response to Physiological and Nonphysiological Concentrations of Insulin.

We employed insulin treatment regimens that mimic normo- and hyper-insulinemic states to mirror in vivo conditions as nearly as possible. Specifically, for capturing phosphorylation dynamics under normo-insulinemic states, we stimulated cells with 0.1 and 1.0 nM insulin, which are usually associated with normal fed-fast cycles (22). Furthermore, we wanted to assess the kinetics of signaling through kinases, which are downstream and govern effector functions of both the metabolic and growth-factor arms, and hence are considered nodal in the cascade (Fig. 1).

Fig. 1.

Iterative experimental–mathematical approach reveals distinct IS kinetics. (A) Schematic of the insulin/IGF signaling pathway. Components involved in metabolic and mitogenic arms are shown in blue and green, respectively. Phosphorylations measured in this study are highlighted in red. (B) Experimental paradigm and workflow for assaying signaling in response to one-step stimulation. (C) Representative blots for levels of pAKTT308, pAKTS473, and pERKT202/Y204 following insulin stimulation, as indicated. Total AKT/ERK and actin were used for normalization (). (D–F) Quantitation for temporal changes in pAKTS473 (D), pAKTT308 (E), and pERKT202/Y204 (F) from experimental data shown in C. (G) Schematic for simulation model for IS pathway. Components colored green and pink are considered in the model. Refer to for ODEs corresponding to individual reactions (re 1 through 27) shown. (H–J) Quantitation for temporal changes in pAKTS473 (H), pAKTT308 (I), and pERKT202/Y204 (J) from mathematical simulations using differential equations. (K) Quantitation for temporal changes in phosphorylations from mathematical simulations obtained without considering experimental data for pGSK3βS9. Fits only shown for the time points indicated. Value for pS6KT389 scaled down by a factor of 100 for visualization. (L) Prediction of pGSK3βS9 kinetics from the parameters used to compute K. (M) Kinetic behavior of phosphorylated pAKTT308 and pAKTS473 molecules at 1 nM insulin from stochastic simulations. The band represents SD. Fold changes in D through F for each insulin concentration are with respect to their own 0-min time point. All data presented are mean ± SEM (n = 4, n = 4).

Given the importance of liver in modulation of insulin action and integration of whole-body physiology, we employed primary hepatocytes isolated from mice livers. In order to mimic a fed state, we assessed signaling in response to both glucose and insulin inputs as described in Fig. 1. Treating cells with low- or high-glucose–containing media alone, in the absence of insulin, did not elicit a significant response and therefore allowed us to assess nutrient-independent but insulin-dependent changes in signaling (). Additionally, our paradigm ensured that the kinetic evaluation did not have any bearing from residual signals (), and we assessed for phosphorylations until 120 min given the normal dynamics of IS in vivo (31). Moreover, such analyses could be confounded by differential specificities and affinities of antibodies, and variations in the unstimulated phosphorylation signal intensities across treatment regimes. We negated these by comparing and normalizing the phosphorylation signals at time zero across multiple treatments and experimental trials (). Furthermore, changes in the extent of phosphorylation was computed as a fold change across regimes for the same antibodies and not across antibodies.

One-nanomolar insulin treatment led to a rapid activation of downstream AKT signaling, as anticipated for a fed response (Fig. 1 and ) and reported earlier (3, 21, 32, 33). Expectedly, overall signal intensities (i.e., area under the curve) for all the phosphorylation events scored in our assay were positively correlated with increasing insulin concentration from 0.1 to 100 nM (). We define true discovery rate as a statistical measure to compare changes in phosphorylations across time and insulin concentrations. This parameter further validated the robustness of our measurements ().

We found that temporal changes in phosphorylation of nodal kinases, AKT and ERK, were markedly different across insulin treatment regimens (Fig. 1 and ), which was consistent with similar observations by Humphrey et al. (14). This posited nonlinear signal flow through metabolic and growth-factor arms. Such nonlinear and nonmonotonic association of signal intensities w.r.t insulin concentrations was seen for components across the cascade, both in terms of extent of phosphorylation and temporal behavior. For example, activation–inactivation kinetics was starkly different for AKT (T308 and S473) and ERK. While the amplitude of AKT phosphorylation decayed by 120 min, ERK phosphorylation showed a distinct second wave of activation (Fig. 1 and ). Similarly, initial induction of phosphorylation of GSK3β and S6K was phase-delayed in response to 0.1 and 1.0 nM insulin treatments and continued to remain elevated long after phosphorylation on AKT started to extinguish (). While delayed response by GSK3β and S6K is consistent with their downstream functions, our paradigms allowed us to delineate signaling kinetics under both normo- and hyper-insulinemic conditions.

Furthermore, we found that maximal phosphorylation and its sustenance varied drastically for the nodal kinases AKT and ERK between normo- and hyper-insulinemic paradigms. While supraphysiological concentrations of 10 and 100 nM insulin induced a rapid, high amplitude peak for AKT phosphorylations at both T308 and S473, the signals seemed to extinguish faster (Fig. 1 ). On the contrary, in response to normo-insulinemic treatments, signal sustenance was longer despite lower peak amplitude. ERK, on the other hand, displayed a flip vis-à-vis very-rapid increase in phosphorylation under hyper-insulinemic states, which was not observed in response to 0.1 and 1 nM insulin treatments (wherein there was a decrease in pERK at very early time points) (Fig. 1). These indicated insulin concentration–dependent programming of signaling kinetics and prompted us to investigate the parameters that could define such behaviors. Therefore, we set out to elucidate kinetic determinants of signal architecture using both experimental and mathematical approaches.

Mathematical Modeling of Signaling Kinetics and Predictive Assessment of Key Phosphorylation–Dephosphorylation Dynamics.

We mathematically modeled the signaling cascade using our experimental results with an aim to extract kinetic parameters that define the network. The details of model construction, deduced parameters and their validation are described in . Briefly, the model was constructed based on the rationale of using minimum-free parameters (Fig. 1). Parameters for simulation were obtained by comparing the experimentally measured kinetics of protein components (termed as state variables): pAKTT308, pAKTS473, pS6KT389, pGSK3βS9, and pERKT202/Y204 at time points 0, 5, 15, 30, 60, and 120 (min) (). Robustness of the simulation trajectories was checked taking into account variables such as initial protein concentrations and loading deviations to parameters across time points. As depicted in ( and detailed in ), neither varying the initial protein concentration nor addition of loading deviations to the parameters significantly alter the simulated phosphorylation dynamics. We considered the IS network as a set of coupled biochemical reactions and used ordinary differential equations (ODEs) to describe the system (respective components and biochemical reactions of the model are detailed in ). Optimized and fitted values (reaction rate constants and initial concentrations of protein components) and predicted or derived variables (dynamic range, peak intensities and noise) are together termed as “kinetic parameters.”

As shown in Fig. 1 and , the simulation results for pAKTT308, pAKTS473, pS6KT389, pGSK3βS9, and pERKT202/Y204 were consistent and nearly overlapping with the experimental data across insulin concentrations. Notably, the differential between experimentally and mathematically derived values was smaller or comparable to the inherent deviation among experimental sets for most components assessed across time points and insulin concentrations (). Furthermore, the magnitude of error in simulations was also found to be less than the experimental error (). Taken together, these showed a high degree of concordance between the experimental and simulated phosphorylation dynamics.

To evaluate effectiveness of our parameters and their dependence on insulin concentration, we simulated phosphorylation dynamics by altering the input insulin concentrations but without changing the parameters for any given concentration. As can be seen in , changing [Ins] (reaction 1 and 2, as per ) led to significant deviations from experimental data.

Next, we simulated the entire system of reactions without considering the experimental data for pGSK3βS9 as one of the ways to check our mathematical model (detailed in ). This also now allowed us to independently predict the kinetics of pGSK3βS9. We found that most of the parameters were comparable (), and the simulated dynamics of all the five components were consistently similar to the ones that were obtained earlier (Fig. 1 and ).

Additionally, we also used stochastic simulations to provide an alternative approach to validate our mathematical predictions, which qualitatively resembled the deterministic approach for quantities as described in Fig. 1. Apart from predicting the response at a population level, this allowed us to determine the fluctuations in the system and compare it with signaling topology or cascade architecture (see Fig. 3 ).

Fig. 3.

Kinetic gates and maximum connectedness are associated with robust topology under normo-insulinemic states. (A) KON/KOFF ratios for phosphorylations across the signaling cascade. Insulin concentrations of 0.1 to 10 nM are depicted separately. Events corresponding to numbers on x-axis are depicted in the table along with the respective reactions they correspond to. Differential equations corresponding to individual events are detailed in . The yellow band refers to the kinetic barrier applied between 0.1 and 10. (B) Noise in signal for phosphorylation at pAKTT308 and pAKTS473 in response to different insulin concentrations, as indicated. (C) Correlation matrix depicts degree of relatedness between phosphorylation events and their evolution with increasing insulin concentration. (D) Network analysis depicting degree of connectedness across insulin concentrations, as indicated. Dashed line represents negative correlation. Significance in correlation: White (P < 0.05) < Blue (P < 0.005) < Green (P < 0.0005) < Yellow (P < 0.00001) as observed by Student’s t test. (E) Number of edges and nodes in a 1 nM network substituted with 10 nM values, as indicated. (F) Network maps of 1 nM insulin perturbed with 10 nM pAKTT308 (a) and pAKTS473 (b).

Furthermore, using optimized parameters, obtained by comparing simulations with experiments as input, we mathematically deduced the activation dynamics of mTORC1 and mTORC2 complexes (). Even though TOR is the common kinase active component, mTORC1 and mTORC2 complexes seem to display distinct properties in terms of upstream activatory cues, substrate specificity, and downstream effects (7, 8). Intriguingly, kinetics of mTORC1 and mTORC2 were starkly different w.r.t each other and independent of the kinetics of downstream phosphorylation events (pS6KT389 and pAKTS473) ( and Fig. 1).

Such analysis also allowed us to assess concordance (or lack of it) between upstream primary kinase and downstream phosphorylation. In the context of mTORC2, which is the primary kinase for AKTS473 phosphorylation, we found discordant phosphorylation dynamics between the two components, which has been observed earlier (34–37). This highlights the presence of additional regulatory steps in controlling activation/inactivation of AKT (Fig. 1 and ). We found similar disparate responses between mTORC1-kinase and its substrate S6K () and AKT-GSK3β (described in Fig. 2), which were unexpected. These clearly indicated that upstream and downstream phosphorylation dynamics are decoupled and hinted at the existence of kinetic gates (see Fig. 3).

Fig. 2.

Continuous variable parameters and high-/low-pass filters determine concentration–dependent insulin response. (A) Extent of change in phosphorylation at pAKTT308 and pAKTS473 across time points between 0.1 and 1 nM insulin. (B) Phase diagram depicting relationship between peak intensity and decay time from simulated data. For a few points, the error bars are smaller than the size of the symbols. (C) Estimated peak amplitude for signaling components as a function of varying insulin concentration depicts dynamic range. Values corresponding to pS6KT389 are scaled down by a factor of 100 for visualization. (D) Dynamic range for signaling components computed from C, as labeled. All data presented are mean ± SEM (n = 4). P value is as observed by Student’s t test.

AKT-Dependent Responsivity to Insulin Is Determined by Phosphorylation at S473.

Several reports have highlighted the necessity of dual phosphorylation of AKT at T308 and S473 for maximal activity (32, 38, 39). Despite this, their relative contributions in terms of gain in signal strength and response to dynamic changes in insulin concentrations are still unclear. Thus, we computed percentage gain in signal against physiological insulin concentrations of 0.1 and 1 nM for pAKTT308 and pAKTS473 (Fig. 2). Comparatively, dose responsiveness and temporal dynamics of pAKTS473 were more pronounced than pAKTT308 (Fig. 2 ). However, it was also intriguing to find that pAKTT308 showed a markedly high response to hyper-insulinemic (10 nM) inputs (Fig. 2). Together, these illustrate differential regulation of phosphorylations at AKTS473 and AKTT308 under normo- and hyper-insulinemic states.

Normo-insulinemic Inputs Follow Nonconcordant Activation and Decay Kinetics of Nodal Signaling Events.

Threshold of activation and dynamic range are key determinants of responsivity in signaling, especially when the inputs are dynamic, as in the case of IS. Hence, we wanted to determine 1) the relationship between peak amplitude and decay kinetics and 2) dynamic range and threshold activation, which collectively dictate the physiological output. A phase diagram depicting peak amplitudes and decay time of phosphorylation events highlighted nonconcordance between these for AKT but not for downstream kinases GSK3β and S6K (Fig. 2).

Most experimental approaches in the past have assayed for signaling in response to very-high insulin inputs, which is rarely physiological. Such deterministic evaluation of signaling also masks threshold kinetics, which is critical to encode biological response. Toward this, we mathematically determined phosphorylation dynamics (including at peak and at 120 min post insulin stimulation), across concentrations from 0.1 nM to 1.0 nM, using three independent simulation approaches (Fig. 2 and ). With regards to peak intensities, which contribute to threshold activation, we found parametric scaling matched best with our experimentally determined levels of phosphorylation for not only 0.1 and 1 nM insulin but also for intermediate concentrations of insulin tested (). Notably, we found disparate dynamic ranges for peak activation (Fig. 2 ). For example, despite showing a difference in terms of extent of phosphorylation (peak intensities) (Fig. 2 ), pAKTS473 and pAKTT308 displayed large and nearly overlapping dynamic range (Fig. 2). On the other hand, pS6KT389 and pGSK3βS9 reached saturation at lower concentrations of insulin (Fig. 2 ).

Since our simulations predicted nonsaturation dynamics for pAKT, we wanted to experimentally verify whether this was indeed the case. We specifically chose 0.3 nM and 0.6 nM, as the response is linear at 0.3 nM and begins to plateau at 0.6 nM insulin. As shown in , our experimental results were consistent with the mathematical predictions and clearly indicated that pAKTS473 indeed displayed a large dynamic range to insulin inputs. Furthermore, EC50 for pAKTS473 was 0.82 nM of insulin () and consistent with earlier reports (17). Now, we show that saturation dynamics (based on peak amplitude) are distinct from EC50 (based on the total signal intensity or area under the curve), highlighting the importance of temporal behavior of both the extent of phosphorylation and signal persistence. Taken together, nonconcordance between peak activation and decay kinetics across signaling components raised the possibility of existence of 1) kinetic insulation of signals and 2) memory of fasted-insulin inputs, which together would define the fed-insulin response.

Diverse Insulin Inputs Generate Differential Kinetic Gates and Signal Noise.

In a multicomponent and multistep signaling cascade, as in IS, it is important to determine parameters that 1) exert control over information flow through the network and thus define topology and 2) those that maintain robustness of the network itself. Kinetic insulation or gating, a selective temporal bias of phosphorylation, has been proposed as one of the key determinants of nonuniform flow of information as inferred by mathematical approaches on synthetic cellular signaling cascades (1, 30). However, it has not been applied to complex and dynamic signaling systems such as IS under physiologically relevant states.

In this context, we used our experimental data and mathematical simulations () to deduce kinetic gates that define architecture of IS. Simplistically, biased retention of specific phosphorylations could be achieved by altered rates of appearance (KON) and decay (KOFF). Thus, to reveal kinetic gating/insulation, we computed rate constants for phosphorylation events, which included known feed-forward and feed-back regulatory inputs (Fig. 3 and ). Furthermore, we assume that very-high or -low ratios of KON/KOFF would constitute kinetic “gates.” We identify a “gate” in a reaction if the KON/KOFF ratio is very high or very low beyond a factor of 10. This corresponds to the situation in which the effective free-energy/-chemical potential is more than twice compared to the thermal energy (10 ∼ e2.3).

KON/KOFF ratios were determined for key phosphorylation events in IS, labeled as events 1 through 12 to capture phosphorylation dynamics of key components and corresponding to respective phosphorylation/de-phosphorylation reactions 1 through 27 as described in Fig. 3 and . Plotting KON/KOFF ratios clearly indicated distinct kinetic barriers across insulin concentrations (Fig. 3). The existence of input-dependent kinetic barriers was corroborated by insulin concentration dependence of parameters that govern phosphorylation dynamics ().

Notably, at insulin concentration of 1 nM, which mimics a physiologically fed state, most reactions were not gated and were unlike the response to very-low and very-high insulin concentrations. We also observed a peculiar pattern between upstream priming events {viz. events 1 [Ins+IR ⇄ p1IRC (phosphorylated insulin receptor complex 1)], 6 [ppAKT+mTORC1 ⇄ ppAKT+pmTORC1], 10 [ppAKT+GSK3β ⇄ ppAKT+pGSK3β], and 11 [p1IRC+Raf ⇄ p1IRC+Raf*]} and their effector or downstream phosphorylations {events 4 [p1IRC+pAKTS473 ⇄ p1IRC+ppAKT], 5 [pmTORC2+pAKTT308 ⇄ pmTORC2+ppAKT], and 12 [Raf*+ERK ⇄ Raf*+ppERK]} vis-à-vis kinetic gates specifically at 0.1 nM (Fig. 3). While the reactions for upstream priming events were associated with a positive gate, their downstream targets have negative gates. Although speculative, such low thresholds of activation at fasted-insulin concentrations of 0.1 nM can possibly ensure high sensitivity to limiting insulin doses while a higher KOFF for effector kinases prevents a run-away response. Conversely, we observed that there was a distinct uncoupling between the metabolic and growth-factor arms at supraphysiological concentrations of 10 nM with hyper-activation of cues for macromolecular synthesis (S6K activation and GSK3β inhibition, events 8 and 10, respectively), growth, and proliferation (pERK formation, event 12).

Noise in biology is often considered as being important for mounting a robust response in addition to generating functional heterogeneity and flexibility (33, 40, 41). Since phosphorylation of AKT is one of the central events that is used as surrogate for IS, and given the differential dynamics of pAKTT308 and pAKTS473, we wanted to assess their individual contributions to functional flexibility. Therefore, we computed noise in their signaling (fluctuations on mathematically determined concentrations; see ). As shown in Fig. 3, we observed that lower concentrations of insulin generate more noise than higher concentrations across time points assessed. This could be attributed to a log normal distribution (for early time points) and stationary distribution for later time points, when mean is saturated (). While being in general agreement with similar measurements of other biological parameters, this also indicated that the differential phosphorylation dynamics of AKTT308 and AKTS473 were independent of noise. However, it was also surprising to find that for lower insulin concentration of 0.1 nM, noise associated with both pAKTT308 and pAKTS473 decreased over time. Given that the activation profile of AKT at 0.1 nM is slower and takes longer to peak, we speculate that this reduction in noise is possibly due to stabilization of activated pAKT molecules at later time points. Conversely, early peaking of signal for higher insulin concentrations could stabilize the molecules much earlier to cause an overall dampening of noise.

Robust IS Topology Is Achieved at Physiological Insulin Inputs.

Degree of connectedness between components within a network is used as one of the measures to assess topology and robustness, and is defined by how correlated their responses are as a continuum. Therefore, we set out to ask whether supra/physiological inputs of insulin had any bearing on signaling topology.

Considering the nonmonotonic and nonlinear nature of phosphorylation dynamics, we employed multiple tests to compute correlation matrices for different insulin concentrations (Fig. 3 and ). Despite potential caveats associated with such correlations for a network that changes over time, our findings do point out a better concordance between the entities measured under physiologically relevant settings.

Next, we checked whether a high degree of correlation in response to fed- (1 nM) and fasted- (0.1 nM) insulin inputs had any impact on connectedness within the network. For a maximally connected network of “n” nodes, the maximum number of edges would be n(n-1)/2, while the minimum number of edges would be (n-1). Applying this to a five-component system (as in our case) gives a theoretical upper bound with 10 connections, although 34 undirected nonisomorphic graphs can be realized. We found that when five nodes corresponding to quantities pAKTS473, pAKTT308, pS6KT389, pGSK3βS9, and pERKT202/Y204 were used, maximum connectivity was obtained at physiological concentrations of insulin (at 0.1 nM and 1.0 nM) (Fig. 3 and ).

Distinctly, the degree of connectedness decreased at 10 nM insulin and the node corresponding to pERK was disconnected. More interestingly, this was associated with a flip in pERK dynamics between 0.1, 1, and 10 nM insulin concentrations (Fig. 1). In this context, deducing change in flux of pERK (in comparison with pAKTT308 and S473) revealed a dramatic shift indicative of altered phosphorylation/de-phosphorylation dynamics in response to supraphysiological inputs of insulin (). Together, these posited a decoupling of the metabolic and growth-factor arms, which could have a bearing on pathophysiological implications associated with hyper-insulinemia.

This prompted us to check which of the measured quantities from hyper-insulinemic regimes perturb the degree of connectedness in the network and impinge on coupled information flow under pathological states. Toward this, we substituted values for each of the nodes from 10 nM measurements in the 0.1 and 1 nM insulin networks. Perturbation of every component led to a reduction in both the number of nodes as well as edges (Fig. 3 and ). Notably, while perturbation of pAKTS473 caused disappearance of some edges, perturbation of pAKTT308 brought down the connections drastically from 6 to 2 (Fig. 3 ), indicative of a broken network.

Pulsatile-Fasting Insulin Rewires Response to Fed-Insulin Inputs Akin to Memory.

Uniquely, insulin is released in a pulsatile manner during a fasted state (42, 43), which is followed by a biphasic secretion in response to fed nutrient inputs. As mentioned earlier, while most studies on signaling dynamics have used high concentrations of insulin, there are no reports that have investigated kinetics and topology vis-à-vis pulsatile fasted-insulin inputs. Moreover, if/how a low-dose, fasted-insulin input shapes signaling architecture in response to a fed-state bolus insulin secretion has not been addressed thus far.

To this end, we pulsed hepatocytes with 0.1 nM insulin and then subsequently treated with 1 nM insulin as a proxy to physiological dynamics of fasted- and fed-insulin inputs, as indicated (Fig. 4). Surprisingly, we found that there was neither sustenance nor an enhanced response to consequent insulin pulses for both pAKTT308 and pAKTS473 (Fig. 4 and ), which was unanticipated. This striking loss of pAKT signal by the end of fourth pulse (at 0’) was distinct from a continuous step treatment as described earlier (Fig. 1 ) and indicated a memory of signaling. This behavior was not seen for pERK (Fig. 4 and ). pAKT levels also reached a new baseline following pulsatile insulin stimulation (Fig. 4 ). This new reset point of pAKT also changed the kinetics following 1 nM insulin treatment, which was distinct from ERK, GSK3β, and S6K phosphorylations (Fig. 4 and ). Pulsing with Earle’s Balanced Salt Solution (EBSS) or low-glucose–containing medium, in the absence of insulin, acted as controls and allowed us to assess signaling changes contributed by insulin (). These results clearly indicated that fasted-insulin pulses created a memory to possibly enhance the response to fed-insulin inputs.

Fig. 4.

Pulsatile-fasting insulin rewires response to fed-insulin inputs akin to memory. (A) Experimental paradigm for mimicking pulsatile-fasted (0.1 nM) and bolus-fed (1 nM) insulin stimulation. P1 through P4 indicate 0.1 nM insulin pulses. (B–D) Quantitation for temporal changes in phosphorylations at pAKTT308 (B), pAKTS473 (C), and pERKT202/Y204 (D) following insulin treatment as in A) (E) pFOXO1S256 dynamics in response to 1 nM insulin treatment, as predicted by mathematical simulation (blue) and subsequent experimental validation (green; refer to Western blots in ). (F) Representative images for immunofluorescence staining of FOXO1 localization following a single-step stimulation with 1 nM or following adaption with 0.1 nM insulin pulses as shown in A. (Scale bar, 50 µm.) Representative cells and respective nuclei are marked with dotted line. (G) Quantitation of FOXO1 from (F) (n = 2, n = 25 to 30 cells for each time point). Inset in G depicts time taken for Nuclear/Cytoplasmic ratio to become less than 50% of 0-min time point. P value is as observed by Student’s t test. (H) Heat maps for changes in gene expression downstream to IS in response to constant 0.1 nM and 1 nM and 1 nM insulin following fasted-insulin inputs (1 nM adapted, as in A) (n = 2, n = 3). Asterisk depicts P values (*P < 0.05, **P < 0.005, and ***P < 0.0005) as observed by Student’s t test. Global color scale depicts fold change with respect to unstimulated cells at 0-min time point. (I) Network analysis showing connectivity among signaling components after treatment with 1 nM insulin following fasted-insulin inputs, as in A. Significance in correlation is White (P < 0.05) < Blue (P < 0.005) < Green (P < 0.0005) < Yellow (P < 0.00001) as observed by Student’s t test. (J) KON/KOFF ratios for phosphorylations across the signaling cascade post 0.1 nM pulsatile insulin adaptation. Numbers on the x-axis represent phosphorylation events as detailed in Fig. 3 and . The yellow band refers to the kinetic barrier applied between 0.1 and 1.0. Fold changes in B through E for each concentration are with respect to their own 0-min time point. Data presented are mean ± SEM (n = 4, n = 4).

IS encodes a gene expression program via downstream transcription factors with FOXO1 acting as one of the key mediators (44). In this context, we used our simulation–experimental approach to predict the kinetics of this downstream-most component and investigate the impact of repeat/pulsatile insulin inputs on its subcellular localization. Simulated FOXO1 dynamics by including two additional phosphorylation–dephosphorylation reactions in our existing model did not change the parameters significantly compared to our original predictions (). Furthermore, these simulated/predicted AKT-dependent phosphorylation dynamics of FOXO1 and subsequent experimentally determined pFOXO1S256 kinetics were nearly overlapping (Fig. 4). Since AKT-dependent nuclear exclusion of FOXO1 is well-established in the field (45), we investigated whether kinetic memory of pulsatile stimulation impacted FOXO1 localization. Change in FOXO1 localization in response to 1 nM insulin was significantly different in a pulsed adapted system, wherein nuclear to cytoplasmic relocalization was observed as early as 5 min when compared to 25 min in a nonpulse adapted system (Fig. 4 and ). Consistent with the phosphorylation dynamics, pulsing in the absence of insulin did not have any bearing on FOXO1 localization (). This clearly indicates that a memory of 0.1 nM pulsatile insulin stimulation either potentiates or sensitizes the kinetics of FOXO1 nuclear exclusion.

Next, we assayed for the transcription of genes downstream to IS (46). As shown in Fig. 4, consistent with the memory hypothesis, 0.1 nM insulin pulses reprogrammed the gene expression in response to 1 nM insulin stimulation.

These prompted us to check whether memory was associated with a higher degree of connectedness and/or gating of signals in the network. We indeed found that the network was more robust following pulsatile-adapted 1 nM insulin stimulation (as compared to 1 nM alone) (Fig. 4 and ). Furthermore, in comparison to 1 nM insulin alone (Fig. 3), pulse priming with 0.1 nM insulin caused selective gating of signals with significant changes in KON/KOFF ratio for pAKTS473 (Fig. 4, reactions 3 and 5).

Repeated Stimulation by Fed Insulin Abrogates the Synergy between the Metabolic and Mitogenic Arms of Signaling.

Although continuous exposure to higher levels of circulating insulin is known to cause resistance and thus metabolic diseases, the kinetic basis for such signaling has not been investigated. In this context, we repeat stimulated hepatocytes with 1 nM insulin, as indicated in Fig. 5. This led to an anomalous response vis-à-vis both metabolic (AKT) and mitogenic (ERK) arms of signaling. While pAKT levels decreased drastically, peak amplitudes of pERK increased following repeated stimulation (RS1 and 2) of fed-insulin inputs (Fig. 5 and ).

Fig. 5.

Repeated stimulation by fed insulin abrogates the synergy between the metabolic and mitogenic arms of signaling. (A) Experimental paradigm for repeated stimulation with 1 nM insulin. (B–D) Quantitation for temporal changes in phosphorylations at pAKTT308 (B), pAKTS473 (C), and pERKT202/Y204 (D) following repeated insulin stimulation. Fold changes for each concentration are with respective to their own 0-min time point. Data presented are mean ± SEM (n = 4). (E) Representative images for immunofluorescence staining of FOXO1 localization following a repeated stimulation with 1 nM insulin as shown in A. (Scale bar, 50 µm.) Representative cells and respective nuclei are marked with dotted line. (F) Quantitation of FOXO1 from (E) (n = 2, n = 25 to 30 cells for each time point). (G) Heat maps for changes in gene expression downstream to IS in response to repeated insulin stimulation (n = 2, n = 3). Asterisk depicts P values (*P < 0.05, **P < 0.005, and ***P < 0.0005) as observed by Student’s t test. Global color scale depicts fold change with respect to unstimulated cells at zero-minute time point. (H) Network analysis showing connectivity among signaling components after repeated insulin stimulation, as in (A). Significance in correlation is White (P < 0.05) < Blue (P < 0.005) < Green (P < 0.0005) < Yellow (P < 0.00001) as observed by Student’s t test.

Aberrant nuclear retention of FOXO1 is associated with insulin resistance (47). As expected, we found a rapid cytoplasmic export of FOXO1 following a single-step stimulation with 1 nM insulin (Fig. 5 and ). However, consistent with a failure to induce AKT phosphorylation following 1 nM repeat stimulation (RS1 and RS2) (Fig. 5 ), FOXO1 was retained in the nucleus (Fig. 5 and ). This abnormal retention was comparable with its localization that was observed following decay in IS in our single-step stimulation paradigm and together indicated resistance in the repeat stimulation paradigm. Moreover, we also observed that repeated stimulation of 1 nM insulin led to altered gene expression profiles, which seemed uncorrelated vis-à-vis genes associated with cellular growth and metabolism (Fig. 5).

Network analysis following this paradigm showed complete loss of connections among signaling components (Fig. 5 and ). This was also apparent with the dynamics of pGSK3β and pS6K, which remain up-regulated despite a down-regulation in AKT signaling ().

Discussion

Coupling nutrient inputs to cellular metabolism, survival, and growth is intrinsically dependent upon IS. Hypo- and hyper-activation of IS leads to various patho-physiologies including diabetes, accelerated aging, and cancer, which are attributed to under- or overphosphorylation of certain IS components (10–13). In this study, using mathematical and experimental approaches, we provide fundamental insights into kinetic parameters that dictate emergent properties of IS, govern its architecture, and couple the metabolic and mitogenic arms under various physiological contexts.

Even though recent phospho-proteomic analyses have unraveled the temporal changes (17, 18), the extent to which phosphorylation dynamics encode information as a function of insulin concentration and/or time is still unclear. We have elucidated the dynamic range of key phosphorylation events under physiologically relevant insulin concentrations (0.1 to 1 nM) as in during fed-fast cycles. This is significant, since in addition to net gain in specific phosphorylations and EC50 (17), dynamic range has been proposed to be a better predictor of cellular response.

Our results based on phosphorylation dynamics at 0.1 and 1 nM insulin also raised the possibility of pAKTT308 and pAKTS473 acting as low-pass and high-pass filters with one of them being a permissive cue, despite having overlapping dynamic ranges. Specifically, while the gain in peak response and kinetics of pT308 was independent of input strength, that of pS473 correlated well with change in ligand concentration under normo-insulinemic ranges but not for hyper-insulinemic concentrations. These results also raise the possibility that while phosphorylation at pT308 may serve to prime AKT, pS473 may determine the extent of overall activation, which could be elicited differentially across insulin concentrations. It is important to note that although Fujita et al. (48) had proposed AKT as a low-pass filter, we now delineate the differential kinetic contributions of phosphorylations at T308 and S473.

Noise in biology is generally regarded to be beneficial for regulating functional flexibility and has been well studied in the context of gene transcription. Given limited knowledge in this regard for signaling cascades (especially for IS), we checked for input versus variance in signal response for the nodal kinase AKT. It was interesting to note that noise in signaling was apparent at physiological concentrations of insulin (0.1 to 1 nM) while it was substantially diminished in hyper-insulinemic regimes. This hinted toward reduced plasticity in signaling under hyper-insulinemic states.

Signal stratification is crucial for sustenance of downstream information even upon input extinction. We discovered that signals are stratified with differential gating in an insulin concentration–dependent manner. These bring to the fore the need to address mechanisms that contribute to these kinetic barriers by affecting KON/KOFF ratios of phosphorylation events, in the future. Emergence of both positive and negative barriers, which are sensitive to changes in flux, posits respective de-/phosphorylation events as key regulatory nodes, especially in response to hyper-insulinemic inputs and during 0.1 to 1 nM insulin transitions (pulsatile regime). Furthermore, we propose that such components would be very attractive candidates for therapeutic interventions to regulate IS and maintain network properties differentially.

Besides kinetic gating, connectedness between signaling components is one of the key determinants of network topology. Despite several studies on signaling cascades across biological systems, little is known about whether/how these parameters contribute to topology, except in cases in which simulations have been carried out for artificial signaling systems. Our study illustrates that maximum connectivity between the signaling nodes, which is often used as a measure of network robustness, is achieved at physiological concentrations of insulin. Conversely, the degree of connectedness is low under hyper-insulinemic states. Importantly, we also underscore the significance of each of the phosphorylations in maintaining the robustness of the topology under normo-insulinemic states. For example, our results suggest that pT308 is a key predictor of network topology.

Others and we have found that metabolic cues under fasting conditions elicit anticipatory molecular mechanisms to mount an efficient fed response (49–51). In this context, we have unraveled that pulsatile 0.1 nM insulin rewires signaling dynamics in response to 1 nM inputs. It was striking to find enhanced net gain in phosphorylation of some (pAKTT308 and S473 and pGSK3βS9) but not all components, akin to memory or anticipation. Conversely, insulin resistance is associated with repeated insulin/nutrient inputs and hyper-insulinemia. In this context, our study also describes kinetic changes in IS dynamics, which can be either causal or consequential to reduced sensitivity under these conditions. Notably, we found that repeated stimulation with 1 nM insulin dampened the AKT response while inducing pERK indicating a disbalance between metabolic and mitogenic arms. This is important because over-activation of either metabolic and/or the mitogenic arm has been described in literature as a driver of metabolic diseases and cancer (52–54). Here, we would like to specifically highlight that the signaling cascade is most robust in response to fed-insulin inputs, which is pulse primed by fasting insulin. Our findings posit that repeated and/or high insulin inputs, including in a clinical setting (inappropriate dosing), could lead to perturbed networks with possible pathological manifestations.

While the paradigms used in our study were designed to capture IS dynamics in primary hepatocytes in vitro, given the complex interplay of IS with other signaling cascades/inputs, kinetics in vivo may well be different. Our study also does not take into consideration dynamics of early phosphorylation events and the endocytosis-based regulation of IS. Even though this is a caveat, it is consistent with current efforts in the field, which is limited by reagents/tools to capture extremely rapid or very early phosphorylation/de-phosphorylation reactions in an endogenous context. Despite advances in the field (such as optogenetic tools used in other signaling cascades), investigating upstream most phosphorylation events remains challenging in an endogenous context. Nonetheless, this study provides insights into insulin concentration–dependent nonmonotonic changes in phosphorylation of nodal effector kinases, besides providing a framework to access the kinetic complexity of information flow in signaling cascades, also in vivo.

In conclusion, our results unravel hitherto unknown kinetic constraints that exert control over components of IS. Notably, we illustrate that these kinetic parameters vary starkly during normo- and hyper-insulinemic states. Given that a discordant signal flow between metabolic and growth-factor arms is associated with diseases, our findings provide fundamental insights into factors that govern this coupling. Our study also raises the possibility of impaired biological outputs in the context of therapeutic interventions using insulin, which have been largely guided by glycemic control. We highlight the importance of discovering regulatory parameters and nodes to complete our understanding of signaling cascades under both normal and pathological conditions.

Materials and Methods

Insulin Treatments.

At 24 h post-plating, the hepatocyte medium was changed to 5% fetal bovine serum containing Dulbecco’s modified Eagle medium (DMEM)–high glucose for 11 h. In order to obtain unstimulated phosphorylation signal (0 min), cells were serum-starved and grown in EBSS for 6 h. Cell culture grade insulin (Catalog No. I0516) with stock concentration of 10 mg/mL was serially diluted in EBSS to make working stocks of 0.02 to 20 nM, such that finally, 10 μL diluted insulin was added per 2 mL culture medium (DMEM–low glucose/high glucose) to obtain concentrations of 0.1, 1, 10, and 100 nM, respectively. For one-step insulin stimulation experiments, 0.1 to 100 nM insulin in DMEM–high glucose was added to the hepatocytes, and cells were collected at time points as described in the results. For pulsatile insulin treatments and repeated insulin stimulation, paradigm modifications are mentioned in Figs. 4 and 5. Every media change was preceded with a phosphate-buffered saline wash to remove residual contamination.