ll OPEN ACCESSiScience ArticleSignaling Heterogeneity is Defined by Pathway Architecture and Intercellular Variability in Protein ExpressionDougall Norris, Pengyi Yang, Sung-Young Shin, ..., Lan K. Nguyen, David E. James, James G. Burchfield lan.k.nguyen@monash.edu (L.K.N.) david.james@sydney.edu.au (D.E.J.) james.burchfield@sydney.edu. au (J.G.B.) HIGHLIGHTS Insulin signaling is heterogeneous between cells in the same population The temporal response of signaling components within a cell is highly reproducible Upstream responses (Akt) can only partially predict downstream response (GLUT4) Protein expression variance is a driver of intercellular signaling heterogeneity Norris et al., iScience 24, 102118 February 19, 2021 ª 2021 The Author(s). https://doi.org/10.1016/ j.isci.2021.102118 iScience Article Signaling Heterogeneity is Defined by Pathway Architecture and Intercellular Variability in Protein Expression Dougall Norris,1,2,8,9 Pengyi Yang,1,3,4,9 Sung-Young Shin,5,6,9 Alison L. Kearney,1,2 Hani Jieun Kim,1,3,4 Thomas Geddes,1,3,4 Alistair M. Senior,1,2 Daniel J. Fazakerley,1,8 Lan K. Nguyen,5,6,* David E. James,1,2,7,* and James G. Burchfield1,2,10,* SUMMARY Insulin’s activation of PI3K/Akt signaling, stimulates glucose uptake by enhancing delivery of GLUT4 to the cell surface. Here we examined the origins of intercel- lular heterogeneity in insulin signaling. Akt activation alone accounted for ~25% of the variance in GLUT4, indicating that additional sources of variance exist. The Akt and GLUT4 responses were highly reproducible within the same cell, suggesting the variance is between cells (extrinsic) and not within cells (intrinsic). Generalized mechanistic models (supported by experimental observa- tions) demonstrated that the correlation between the steady-state levels of two measured signaling processes decreaseswith increasing distance from each other and that intercellular variation in protein expression (as an example of extrinsic variance) is sufficient to account for the variance in and between Akt and GLUT4. Thus, the response of a population to insulin signaling is underpinned by considerable single-cell heterogeneity that is largely driven by variance in gene/protein expression between cells. INTRODUCTION Despite our long-held appreciation of heterogeneity at the species level, we have only recently begun to appreciate the significance of intercellular heterogeneity with respect to biology and disease. Recently, studies of single cells have uncovered that cell-to-cell heterogeneity underpins embryonic development (Biddy et al., 2018; Zhao et al., 2019), cellular metabolism (Sharick et al., 2019; Zhang et al., 2018), immunity (Miragaia et al., 2019; Shalek et al., 2013), neuronal function (Chelaru and Dragoi, 2008; Petitpre´ et al., 2018), clonal evolution in cancer (Navin et al., 2011; Nolan-Stevaux et al., 2013; Wang et al., 2014), and multidrug resistance in bacteria and tumors (Balaban et al., 2004; Gorre et al., 2001; Greaves, 2015; Keren et al., 2004; Lee et al., 2014). Such findings have not only furthered our understanding of biological mechanisms but also necessitated a change in the design of therapeutic approaches to accommodate heterogeneous cellular populations in disease. The insulin signaling pathway is central to the regulation of glucose metabolism. Adipose and muscle cells respond to insulin by redistributing GLUT4 from intracellular compartments to the plasma membrane (PM). This occurs in response to the activation of a signaling cascade initiated by insulin binding to its receptor at the PM, which leads to recruitment and activation of class I phosphatidylinositol 3-kinase (PI3K). PI3K- generated phosphatidylinositol-(3,4,5)-trisphosphate (PIP)3 at the PM serves as a docking site for the Ser/Thr kinase Akt, which is activated by phosphorylation at the PM. Akt is a central regulatory node in the insulin signaling pathway, and its activation is necessary and sufficient for insulin-stimulated GLUT4 translocation (Kohn et al., 1996; Ng et al., 2008; Tan et al, 2012b, 2015). Akt-mediated phosphorylation of the RabGAP AS160/TBC1D4 releases vesicles containing GLUT4 from intra-cellular compartments to the PM. This permits the rapid uptake and storage of glucose (Sto¨ckli et al., 2008; Tan et al., 2012a). We have previously used live cell total internal reflection fluorescence microscopy (TIRFM) to unveil hetero- geneity in the insulin-stimulated redistribution of glucose transporter 4 (GLUT4), to the PM in adipocytes (Burchfield et al., 2013). Intriguingly, individual cell responses were non-stochastic and highly reproducible 1Charles Perkins Centre, The University of Sydney, Sydney, NSW 2006, Australia 2School of Life and Environmental Sciences, The University of Sydney, Sydney, NSW 2006, Australia 3School of Mathematics and Statistics, The University of Sydney, Sydney, NSW 2006, Australia 4Computational Systems Biology Group, Children’s Medical Research Institute, Faculty of Medicine and Health, The University of Sydney, Westmead, NSW 2145, Australia 5Department of Biochemistry and Molecular Biology, School of Biomedical Sciences, Monash University, Clayton, VIC 3800, Australia 6Biomedicine Discovery Institute, Monash University, Clayton, VIC 3800, Australia 7Sydney Medical School, The University of Sydney, Sydney, NSW 2006, Australia 8Present address: Metabolic Research Laboratories, Wellcome-Medical Research Council Institute of Metabolic Science, University of Cambridge, Cambridge, CB2 0QQ, UK 9These authors contributed equally 10Lead Contact *Correspondence: lan.k.nguyen@monash.edu (L.K.N.), david.james@sydney.edu.au (D.E.J.), james.burchfield@sydney. edu.au (J.G.B.) https://doi.org/10.1016/j.isci. 2021.102118 ll OPEN ACCESS iScience 24, 102118, February 19, 2021 ª 2021 The Author(s). This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 1 across multiple stimuli, implying that adipocytes possess a predetermined sensitivity to insulin. As impaired insulin-stimulated GLUT4 translocation to the PM is a key feature of insulin resistance (Cline et al., 1999; Fazakerley et al., 2018; Gonzalez et al., 2011; Hoehn et al., 2008; Tremblay et al., 2001), under- standing the basis for this heterogeneity may further our understanding of insulin resistance. Given the central role of PI3K/Akt in GLUT4 trafficking, we hypothesized that intercellular heterogeneity in insulin-stimulated GLUT4 translocation would be a direct consequence of heterogeneity in PI3K/Akt acti- vation. Using a combination of live single-cell imaging, statistical modeling, and generalized mechanistic models we found that (1) distinct aspects of the Akt translocation response to insulin were independently predictive of GLUT4 responses within single cells; (2) this predictive power was weaker than anticipated, suggesting that additional variance is introduced between Akt and GLUT4; (3) the correlation between Akt and GLUT4 responses in the same cell was likely a consequence of the large number of regulatory steps between them; (4) signaling events closer in time and space had amuch greater correlation within cells; and (5) response heterogeneity throughout the population is driven by extrinsic factors (differences between cells) that are stable over hours such as difference in protein expression between cells. These observations are crucial for furthering our understanding of how signaling within single cells contribute to the response of the population and how this relates to the development of disease. RESULTS Single-cell Akt responses to insulin are heterogeneous and reproducible We and others have previously demonstrated that insulin-stimulated GLUT4 translocation to the PM is markedly heterogeneous between individual cells (Burchfield et al., 2013; Lizunov et al., 2015) and repro- ducible within a cell (Burchfield et al., 2013). However, the factors that give rise to this heterogeneity remain unclear. Given that Akt is a central node of the insulin signaling pathway and that the recruitment of Akt to the PM is required for GLUT4 translocation (Ng et al., 2008), we hypothesized that heterogeneity in insulin- stimulated PI3K/Akt signaling could drive the heterogeneity observed in GLUT4 translocation. To test this, we first compared the intercellular heterogeneity observed for insulin-stimulated GLUT4 trans- location and Akt recruitment. Consistent with our previous findings using the pHlourin-GLUT4-tdTomato reporter construct (Burchfield et al., 2013), insulin-stimulated delivery of GLUT4-pHluorin to the PM dis- played a high degree of intercellular heterogeneity (Figure 1A). At the same physiological dose of insulin, Akt2 recruitment was also highly heterogeneous (Figure 1B). Our previous observation that variability in cellular GLUT4 responses were reproducible over repeat stimuli indicates that the variability is non-sto- chastic and represents a predetermined sensitivity to insulin. We next tested whether the heterogeneity in insulin-stimulated Akt2 translocation was also non-stochastic. Cells were exposed to 1 nM insulin for 20 min, followed by a 2-h washout with basal media and a second 1-nM insulin stimulation. The response to the second stimulus was highly correlated with the primary response (Figures 1C and 1D; r = 0.98 versus r = 0.71 for a random pair). Together, these data demonstrate that the recruitment of Akt2 in response to insulin is both heterogeneous and non-random. This is consistent with our observations of GLUT4 traf- ficking (Burchfield et al., 2013) and with the possibility that heterogeneity in insulin-stimulated Akt activa- tion may drive heterogeneity in GLUT4 translocation. Validation of a co-expression system for assessment of insulin-stimulated Akt2 and GLUT4 recruitment To further explore the relationship between Akt and GLUT4 within single cells, we conducted a series of live-cell imaging experiments where Akt2-tagRFP-T andGLUT4-pHluorin were co-expressed in adipocytes. Across five independent experiments in which a total of 206 cells were imaged, we observed extensive het- erogeneity in both the Akt2 and GLUT4 responses (Figure 2A). Individual cells from the five independent experiments were well-mixed in the principal-component analysis (Figure 2B), suggesting the heterogene- ity we observed in Akt2 and GLUT4 responses are unlikely to be driven by experimental batch effects. There was a strong linear correlation between the expression levels (as assessed by epifluorescence) of Akt2- tagRFP-T and GLUT4-pHluorin, indicating that there was no preferential expression of one protein over the other (Figure 2C, upper panel). The implication is that both constructs are delivered at similar levels during co-electroporation. Finally, there was no correlation between Akt2 expression and the GLUT4 response, indicating that the observed GLUT4 responses were not driven by the degree of Akt2 overex- pression (Figure 2C, lower panel). Taken together, these data suggest that the effect of overexpression on Akt2 and GLUT4 responses was minimal. ll OPEN ACCESS 2 iScience 24, 102118, February 19, 2021 Maximum Akt2 recruitment is a predictor of the GLUT4 response in the same cell We predicted that we would observe a high degree of correlation between the Akt and GLUT4 response within a single cell. To explore this, we assessed the major response parameters between Akt and GLUT4. For Akt these were the area under the response curve (AUC), response maximum (Max), and the response profile (Shape; calculated as the AUC of min-max normalized response). For GLUT4 these were the Max, response half-time, and response lag (see Transparent Methods). The strongest correlation we observed was between the Akt AUC and the Max GLUT4 response (r = 0.52, Figures 3A and 3B), followed by the Akt Max and GLUT4 Max (r = 0.47, Figures 3A and 3B). Interestingly, the shape of the Akt response was also predictive of the GLUT4 Max (r = 0.23, Figures 3A and 3B). To further dissect the correlation between Akt Shape and GLUT4 Max, we performed a fuzzy c-means clustering and found two distinct Akt response profiles; one group of cells had a more sustained Akt response after reaching a peak (Cluster 1; Figure 3C), and the other peaked and then sharply declined (Cluster 2; Figure 3C). Cells in cluster 1 had a higher GLUT4 max than those in cluster 2 (p = 2.8x10-3; Fig- ure 3C, bottom panel) suggesting that a more sustained Akt response leads to improved GLUT4 translocation. Figure 1. Akt2 membrane recruitment is heterogeneous and non-random (A and B) The insulin-stimulated (1 nM) translocation of both GLUT4 (A) and Akt2 (B) to the PM is heterogeneous within a population of cells. (C) Akt2 recruitment to the PM is reproducible following a 2-h washout in response to 1 nM insulin within a single cell. (D) Statistical test comparing reproducibility of Akt2 recruitment (in response to 1 nM insulin before and after 2-h washout) within a cell (red boxes) and between cells (light blue boxes). p values were calculated from a two-sided Wilcoxon rank- sum test. The boxes capture the lower and upper quartiles with the median displayed as a horizontal line in the middle; whiskers are 1.5*IQR. ll OPEN ACCESS iScience 24, 102118, February 19, 2021 3 Akt max and shape are independent predictors of the GLUT4 response To further assess the predictive power of various Akt response characteristics toward the GLUT4 maximum we applied a random forest regression model. Consistent with the correlation analysis (Figure 3A), the maximum of the Akt response was highly predictive of the GLUT4 Max, with the Akt response profile (i.e., Akt shape; cluster membership scores) also displaying predictive power (Figure 4A). Interestingly, although the AUC of Akt was highly correlated with GLUT4 Max (Figure 3A), random forest analysis sug- gested that there was very low predictive power of this Akt characteristic toward the GLUT4 response (Fig- ure 4A). As AUC is a combination of both the response magnitude and response profile, we hypothesized that the low predictive power of Akt AUC was due to the high redundancy when these features were also present in the analysis model. To investigate this, we performed partial correlation analyses between GLUT4 max and different combinations of Akt features. Although all three Akt features were strong predic- tors for GLUT4max (Figures 4B and 4C), partial correlation analyses highlighted the redundancy of Akt AUC when Akt max and shape are included (Figure 4B), and, notably, these analyses revealed that Akt max and shape are important independent predictors of GLUT4 max (Figure 4C). Interestingly, the statistical modeling suggests that the combined power of both Akt2 variables (Max, Shape) explains ~27% of the variance in GLUT4 (Figures 3A and 3B, Akt AUC versus GLUT4 max r = 0.52, r2 = 0.27). This suggests that variance is introduced between Akt and GLUT4. Mechanistic modeling demonstrates that intercellular variability in protein expression can account for the increase in variance between Akt and GLUT4 We surmised that variance introduced into the systemwouldmost likely result fromeither intrinsic noise (such as the stochastic nature of biochemical reactions; intracellular variance) or from extrinsic noise (such as the vari- ability in gene/protein expression between cells; intercellular variance). The high degree of reproducibility in Akt and GLUT4 responses within single cells subjected to repeated stimuli (Figures 1C and 1D; Burchfield Figure 2. Quality control analyses of 5 GLUT4 and Akt2 recruitment datasets reveal no deleterious effects of overexpression (A) Temporal GLUT4-pHluorin trafficking and Akt2-tRFPt recruitment in individual 3T3-L1 adipocytes upon 1-nM insulin stimulation visualized by TIRFM. (B) Principal-component analysis plots summarizing response profiles of GLUT4 and Akt2 in each cell into two principle components. Cells are color coded by experiment. (C) Correlations between expression levels of each construct (upper panel) and Akt2 expression level andmaximal GLUT4 response (lower panel). Data are from 206 cells from 5 independent experiments (n = 43, 40, 31, 34, 58 for each experiment, respectively). ll OPEN ACCESS 4 iScience 24, 102118, February 19, 2021 et al., 2013) argues against a significant contribution of intracellular variance driving response heterogeneity. Instead, it suggests the variance is intercellular and is stable within a cell over the timescale we are considering (hours). Given that it is now well established that gene expression in single cells is heterogeneous (Buettner et al., 2015; Evers et al., 2019; Kashima et al., 2018) and protein expression typically correlates with gene expres- sion (Yang et al., 2019) we hypothesized that differences in the expression of individual proteins between cells would introduce stable intercellular variance throughout the pathway. To assess the degree of variance we could expect in the expression of a single gene between cells, we explored the range of gene expression that exists in 18 cell types from the TabulaMuris Consortium (Tabula Muris Consortium, 2018). In these datasets 2,000–5,000 genes were quantified per cell (Figure 5A). Toconservativelyestimate themagnitudethat theexpressionofa singlegenecanvarybybetweencellswecalcu- lated thedifferencebetween themeansof theupperand lowerquartile for eachgeneandcalculated theaverage of each gene for each cell line (Figure 5B).We observed that themajority of genes vary by approximately 32-fold between cells, with a range from1- to 1,000-fold (Figure 5B). For relevance to insulin signal transduction, we then limited the analysis to genes definedby theGO termGO:0007169, transmembrane receptor protein tyrosine ki- nase signaling pathway. Here we observed that expression ranges between 9- and 70-fold with the majority of genes varyingbymore than32-fold (Figure5C). Similar variabilitywasobserved in3T3-L1fibroblasts, theparental cell line used in this study (Figures S1A and S1B). These data demonstrate that gene expression is highly variable between cells and is consistent with protein expression being similarly variable between cells. We next con- structed two kinetic models of a generalizedmulti-tier signaling network (with 6 nodes) that resembles a typical signaling cascade triggered by an extracellular input (e.g. insulin), first without and then with a simple negative feedback (Figures 6A and 6B). The latter model was considered, as several negative feedbacks have been Figure 3. The magnitude and shape of insulin-stimulated Akt2 recruitment are predictive of GLUT4 response heterogeneity in single cells (A) Correlations of intracellular GLUT4 and Akt response characteristics. (B) Scatterplot of individual Akt2 response features with respect to GLUT4 maxima. Fitted lines are from linear regression model. The color points in the third panel denote the clusters of each cell based on (C). (C) Akt recruitment clusters as determined by response shape and their corresponding GLUT4 maxima. For the line plot data are mean G SD. For the boxplot, the boxes capture the lower and upper quartiles with the median displayed as a horizontal line in the middle; whiskers are min and max. See also Figure S4. ll OPEN ACCESS iScience 24, 102118, February 19, 2021 5 described in the insulin signaling pathway (Copps andWhite, 2012; Dibble et al., 2009) and we and others have observedbehaviors consistentwith thepresenceofnegative feedback (Balbiset al., 2000;Cedersundet al., 2008; Ebner et al., 2017; Norris et al., 2017b). A detailedmodel description including model reactions, equations, and parameters is given in the Supplementary Text. Using thesemodels, we performed 50 independent time course simulations at a nominal set of kinetic parameter values, where the expression of each network node (Sn and Pn) was allowed to randomly varywithin a 25-fold range (5-fold upanddownof the respective nominal values), which is in line with the variability in gene expression observed in the TabulaMuris and 3T3-L1 fibroblast datasets (Fig- ures 5A–5C and S1). To avoid bias with the nominal kinetic parameter values, we repeated these simulations for 1,000 randomly sampled sets of kinetic parameters. In general, the linear cascade generated rectangular hyper- bola-shaped response profiles that shifted toward sigmoidal responses for nodes further downstream (Figures 6C, S2A, and S2B). On the other hand, the introduction of a negative feedback resulted in an overshoot behavior in the initial responseprofile (Figures 6E, S2C, and S2D) that transitioned tohyperbolic-shaped responses further downstream (Figures 6E, S2C, and S2D). This resembles the dynamics experimentally observed for Akt and GLUT4 (Figure 2A) following physiologically relevant insulin stimulation (1 nM). As expected, varying the feed- back strength within a parameter set demonstrates that feedback introduces an overshoot, decreases the response magnitude, and increases the response speed (Figure S3). These data show that the models were able to capture a diverse range of behaviors and importantly that variance in protein expression introduced vari- ance into the system at a level similar to what has been observed in biological settings. Irrespective of the presence of feedback, we observed that the variance in the system is compounded at each node, which leads to a substantial increase in the response variability with distance from the primary node (i.e., aS1) (Figure 6I). Consistent with the simulation result, the population responses of Akt and GLUT4 had coefficients of variation of 18% and 57%, respectively, following stimulation with 1 nM insulin (Figure S4). As a consequence of the increased variance, we observed that the correlation of aS1 and down- stream nodes decreased as the number of nodes between them (distance) increased (Figures 6D and 6F). Importantly, this reduction in correlation holds true even when averaged across all 1,000 different random parameter sets (Figure 6H). The amount of variance at a given node is dependent on the overall length of the cascade at the point being considered, whereas the correlation between two nodes is dependent on the number of steps between them as this dictates howmuch variance is accumulated. Consistent with this, when using a within-cascade node, e.g., aS4 as a reference node, model simulations in both models show that correlation with aS4 gradually decreases for nodes that are positioned either further upstream or downstream (Figure 6I). Feedback suppresses variance in a fashion dependent on the feedback strength (Figure 6G), and this is associated with a decrease in the correlation of aS1 with all subsequent nodes (Fig- ure 6G), suggesting that variance may be an important factor in determining the strength of a correlation. Consistent with this, we observed that the level of variance at node 1 was strongly correlated with how aS1 correlated with downstream nodes (Figure S5). Taken together these data suggest that intercellular vari- ability in the expression of individual proteins introduces variability into the system that accumulates across each node but allows for response reproducibility within a single cell. Figure 4. Akt AUC predictiveness of GLUT4 max response can be decomposed into Akt max and Akt shape (A) Predictive power of cluster membership scores from the two fuzzy c-means clusters, AUC, and Akt2 max and Akt response shape. Time points before the insulin stimulation (from 1 to 10) are included as negative controls. Variables are sorted by their predictive power based on a random forest regression model. B, C) Correlation and partial correlation analyses on Akt AUC redundancy in GLUT4 max response prediction (B), and independence of Akt max and Akt shape in GLUT4 max response prediction (C). ll OPEN ACCESS 6 iScience 24, 102118, February 19, 2021 The proximal signaling nodes PDK1 and Akt are highly correlated The mechanistic modeling predicted that events spatiotemporally closer to Akt recruitment would display a greater correlation thanobserved for nodes further away, suchasGLUT4 translocation. To test this,wemeasured the correlation between PDPK1 and Akt recruitment, which both bind PIP3 at the PM upon insulin stimulation (James et al., 1996; Lasserre et al., 2008; Park et al., 2008). PDK1-eGFP and Akt2-tRFPt were co-expressed in 3T3-L1 adipocytes and imaged by TIRFM following stimulation with 1 and 100 nM insulin (Figure 7A). The fold change in insulin-stimulated Akt and PDK1 recruitment was heterogeneous between cells (Figure 7A). Impor- tantly, Akt and PDK1 recruitment were significantly correlated in bothmagnitude and shape at both doses of in- sulin (Figure 7B), withPDK1 recruitment being able topredict between 58%and 81% (r2 = [0.76]2 and [0.9]2) of the Akt response. We next measured the phosphorylation of Akt at T308 by PDK1. Cells expressing Akt-tRFPt were imaged by TIRFM and stimulated with 1 nM insulin. After 2 min the cells were rapidly fixed and immunostained with an anti-phosphoT308 antibody. The same cells were then re-identified and phospho-T308 staining quanti- fied by TIRFM and correlated with Akt recruitment. Despite the low number of cells (due to experimental complexity), the recruitment of Akt was highly correlated with its phosphorylation at T308, with recruitment ex- plainingmore than 75%of T308 phosphorylation (Figures 7C ,and 7D). These correlations are in stark contrast to the27%of theAkt response that canexplain theGLUT4 response (Figure 3) and support themodelingprediction that closer nodes within a signaling network have a greater correlation. DISCUSSION In recent years, our understanding of the importance of cellular heterogeneity and efforts to investigate single- cell biology has soared. Such studies have been highly informative of biological mechanisms that previously eluded discovery using population-averaged approaches (Chi, 2014; Cooper and Bakal, 2017; Wang and Bod- ovitz, 2010), such as the contributions of adaptive and acquired resistance of tumor subpopulations in response to neoadjuvant chemotherapy (Kim et al., 2018). Previous findings by our laboratory have demonstrated that in- sulin-stimulated delivery of GLUT4 to the PMdisplays intercellular heterogeneity (Burchfield et al., 2013), yet the source of this heterogeneity was not clear. In the present study we tested the hypothesis that insulin-stimulated GLUT4 responses would largely be driven by the activity of Akt, a critical node in the insulin signaling pathway. We employed live-cell imaging to determine whether cell-to-cell heterogeneity in GLUT4 trafficking correlated with the indices ofAkt activation: PM recruitment andAkt phosphorylation. Insulin-responsive PM recruitment of Akt2 was highly heterogeneous across a population of adipocytes, yet individual responses were highly repro- ducible following repeat stimuli, consistent with our observations of insulin-stimulated GLUT4 translocation (Burchfield et al., 2013). Themagnitude and kinetics of the Akt2 recruitment response were strong, complemen- tary predictors of the GLUT4 response, consistent with the known importance of the PI3K/Akt node in this Figure 5. Individual genes are expressed at highly variable levels between cells (A) Boxplot of relative frequency of gene expression across 18 cell lines from the Tabula Muris Consortium. (B) Boxplot of the difference in gene expression observed between cells. Boxes capture the interquartile range with the median represented by a horizontal bar and whiskers representing the upper and lower 1.5*IQR values. Points denote outliers. (C) Scatterplot illustrating the expression range of genes defined by GO term transmembrane receptor protein tyrosine kinase signaling pathway. x axis denotes the median gene expression across the 18 cell lines. y axis denotes the median log2 fold change in gene expression. Each point denotes a gene and has been colored by the coefficient of variation (CV). See also Figure S1. ll OPEN ACCESS iScience 24, 102118, February 19, 2021 7 Figure 6. Spatiotemporal correlation between nodes within a generalized multi-tier signaling network (A) A generalized multi-tier signaling cascade with no feedback loop. Reactions (in this case phosphorylation) within each tier are reversible, and the active form of the nodes can be converted to the inactive form by phosphatases (Pi, i = 1 . n). (B) As in (A) but where S1 is under regulation of a negative feedback loop. (C) Representative simulated time course response profiles of each network node in (A) following insulin stimulation at t = 0; each curve represents one of 50 different runs where the expression of network nodes were randomly varied. (D) Scatterplots showing steady-state correlation between the primary upstream node (active S1, aS1) and downstream nodes from the representative simulated time course response profiles from (C). (E and F) (E) As in (C) but for the negative feedback network (B). (F) As in (D) but for (E). (G) Coefficient of variation for each of the nodes across all parameter sets for both models. (H) Correlation coefficients between the simulated steady state of active S1 (primary node) and the downstream nodes showing correlation decreases as their distance increases. The correlation coefficient indicates Pearson’s linear correlation coefficient. The error bar indicates standard deviation from randomly generated kinetic parameter sets (n = 1,000). *p<0.05 compared with aS4/aS4 node with the same feedback, by by two-way ANOVA with Dunnett correction for multiple comparisons. (I) Correlation between the simulated steady state of active S4 and the up- and downstream nodes for both models where the feedback is strong. ****p < 0.0001 for change with distance from node1 by two-way ANOVA. #p < 0.01 compared with no feedback within the same node by two-way ANOVA with Dunnett correction for multiple comparisons. See also Figures S2–S5. ll OPEN ACCESS 8 iScience 24, 102118, February 19, 2021 process. The lack of predictive power of the Akt2 response in the other features of GLUT4 translocation (half- time, lag) suggests that the regulation of steady-state levels of PM-associated GLUT4 is the most important aspect of insulin-stimulated-GLUT4 trafficking, which is not surprising, given this underpins glucose uptake (Klip et al., 2019). Intriguingly, the variability in Akt recruitment could only explain around 25% of the variability in the GLUT4 response. This implies additional sources of variability. Although the variability could either be intracellular (driven by stochastic nature of molecular interactions) or intercellular (driven by differences in gene/protein expression), only an intercellular source of variance can explain the intracellular response repro- ducibility that we observe for Akt and GLUT4. Given that the expression of individual genes can vary by several orders of magnitude between cells, we explored the hypothesis that signaling heterogeneity is driven by vari- ability in protein expression and found that this is likely to be amajor contributor to the variance that is observed for Akt andGLUT4. Themodeling demonstrates that variance in the system is compounded at each node in the pathway. As such, the correlation between nodes decreases as the number of nodes between them increases, which was corroborated by experimental data in our biological system. Together, these data strongly suggest that cell-to-cell difference in protein expression is a major driver of intercellular signaling heterogeneity. As hypothesized, our data suggest that Akt is predictive of the GLUT4 response, consistent with the impor- tance of the PI3K/AKT node in this pathway. Additionally, we observed that variance was introduced Figure 7. Adjacent components of the insulin signaling pathway are highly correlated The magnitude and shape of insulin-stimulated PDK1 and Akt recruitment to the plasmamembrane are highly correlated. (A) 3T3-L1 adipocytes were co-electroporated with Akt2-tRFPt and PDK1-eGFP and treated with 1 and 100 nM insulin. Recruitment was assessed using TIRFM. Each line is the response of an individual cell, colored in accordance with the magnitude of recruitment. (B) Comparisons of the PM recruitment maxima (left) and shape (right) observed for Akt and PDK1 in response to 1 (upper) and 100 (lower) nM insulin. (C) PM recruitment of Akt2-tRFPt in single 3T3-L1 adipocytes following a 2-min 1-nM insulin stimulation, as measured by TIRFM. Each cell’s response was ranked from dark to light blue, in order of increasing T309 phosphorylation quantified by immunofluorescence. (D) Scatterplot indicating the degree of correlation between T309 and Akt2-tRFPt PM recruitment maxima for the 11 cells in (C). ll OPEN ACCESS iScience 24, 102118, February 19, 2021 9 between Akt and GLUT4 that appeared to contribute substantially to GLUT4 heterogeneity. This variance was consistent with it being extrinsic in nature and primarily the result of differences in protein expression between cells. Akt is a critical regulator of insulin-stimulated GLUT4 translocation. Indeed, themost power- ful predictor of the GLUT4 response in our experiments was the area under the curve for Akt translocation to the PM, which could be decomposed into the response maximum and the shape. This suggests that it is both the extent of Akt recruitment and the duration of the Akt response that are important for signaling to GLUT4. The effect of magnitude on a signaling response is quite intuitive, whereas the importance of shape is less so. The 1 nM insulin Akt response profile (shape) broadly separated into two classes, those with a strong overshoot and those with more of a square-wave or sustained response. Overshoot behaviors are consistent with negative feedback (Cheong and Levchenko, 2010). Indeed, when we introduced a negative feedback into our generalized pathway model we observed an overshoot in the upstream responses in many of the randomly sampled parameter sets (Figures 6E, S2C, and S2D), which was not present in the linear cascade. Given that negative feedback exists in the Akt pathway (Copps and White, 2012; Dibble et al., 2009), it is tempting to speculate that the Akt response profile is an indicator of the strength of this feedback and that the reduced GLUT4 response observed in the cells with a stronger overshoot is a result of stronger negative feedback. Irrespective of whether it is driven by feedback, these data suggest that sustained recruitment of Akt to the PM increases signal transmission to its substrates such as TBC1D4/ AS160, which is required for release, tethering, and docking of GLUT4 storage vesicles (GSVs) at the PM (Bai et al., 2007; Brewer et al., 2011; Gonzalez and McGraw, 2006; Xiong et al., 2010). Despite the strong correlation between the Akt response andGLUT4, variance in Akt only accounted for around 25% of variance in the GLUT4Max response (Figure 3) and had an evenmore limited predictive power on other features of the GLUT4 response such as lag and half-time. These data suggest that variance is introduced be- tweenAkt andGLUT4.Consistentwith this, the variance in theGLUT4 responsewas approximately 3-fold higher than for Akt in response to 1 nM insulin (57% vs 18%, respectively). So, what is the source of this variance and where does the variance in Akt itself come from? One possibility is the random nature of biochemical reactions that is driven by the random interaction of biomolecules (intrinsic noise). This sort of stochastic variability would lead to the introduction of noise at each step of the pathway; however, it would not give rise to the high intra- cellular response reproducibility that is observed for Akt (Figure 1D), GLUT4 (Burchfield et al., 2013), and in other systems (Wada et al., 2020). This highlights a preexisting sensitivity within individual adipocytes whereby the source of variance must give rise to heterogeneity between cells and reproducibility within cells (i.e., extrinsic variance). Gene/protein expression is the most commonly cited example of extrinsic variance, and considering that the expression of individual genes can vary markedly between cells within the same population (Figure 5), we propose that this is likely to be the major source of variance in the system we are considering. Based on pre- vious observations in these cells (Su et al., 2019) protein expression would not be expected to change signifi- cantly over the timescale we are considering, leading to reproducibility within a cell. The generalized signaling models demonstrated that variability in protein expression between cells is able to account for the observed heterogeneity in the insulin signaling pathway. The modeling demonstrates that irrespective of the presence of feedback, the strength of the correlation between nodes in a signaling pathway decreases as the distance between them increases. This occurs because the variance accumulates at each node the signal traverses. Importantly, this is also what we observe in the biological system with higher correlations between Akt and PDK1 recruitment/Akt T309 phosphorylation when compared with GLUT4, which is at least several steps downstream of Akt. We also observe increased variability in GLUT4 when compared with the Akt response (Figure S4), which is also consistent with the modeling. We therefore propose that variation in protein expression levels (likely a direct result of gene expression) is a central factor in determining the outcome of a signal transduction in single cells. This combines with the signaling architecture to determine the response output. The loss of information along nodes of intracellular signaling cascades unveiled in our integrative analysis may be compared to the information loss governed by the data processing inequality principle, a key concept in information theory. However, there are distinctions between the two. Although the data pro- cessing inequality principle posits that information is typically lost through the processing of data, often with an aim to simplify it, the reduction of information in our studied system is instead due to accumulation of noise (variance) at each node as the signal traverses. In this sense, the phenomenon can be linked to the law of error propagation that governs how the variance of some population outcome could be derived from the individual variances of the other populations that contribute to it. An implication of this law is that ll OPEN ACCESS 10 iScience 24, 102118, February 19, 2021 variance typically accumulates when there is integration of signal, which is often the case in signaling path- ways. In our studied signaling cascades, the individual populations would be the signaling nodes, with the final variance being the variance of the end (distal) node. However, instead of the distal node integrating signals from multiple nodes, in our cascades signal is sequentially passed from one node to the next, and variance is accumulated at each node, resulting in a high variance at the distal one. In addition, our modeling shows that variance still increases along the cascade even when a negative feed- back is embedded into our network structure (albeit overall at lower levels). The behavior of noise is more predictable for linear cascades, whereas the presence of feedback complicates intuitive reasoning due to non-linearity and signal integration. We anticipate that depending on the specific structure of the added feedback loops, i.e., where the feedback starts and ends, and whether it is positive or negative, the prop- erties of signal and noise propagation in signaling cascade could be very different. In these cases, kinetic modeling and simulation as similarly done here would be highly useful in exploring these emergent behaviors. These observations described herein lead one to question whether the variance in this system serves a purpose. Is there a biological consequence or has biology exploited the existence of this variance? Stochastic resonance is a general phenomenon whereby noise added to a non-linear system improves the performance of the system (McDonnell and Abbott, 2009). Wada et al. report a high degree of intercellular variability in the response of individual myotubes to an electrical stimulus. Using the theory of mutual information, they demonstrated that the intercellular noise carries information and improves the response dynamics of the population (Wada et al., 2020). Consistentwith their observationwe havepreviously demonstrated that adipocytes display a robust dose response to insulin; however, this is not observed at the single-cell level (Burchfield et al., 2013). This is consistent with variance in insulin signaling carrying information in a similar fashion. In our mechanistic models, we observe that the amount of variancewithin node 1 is tightly linked to howwell the node 1 response correlates with downstream nodes (Figure S5), with increased variance leading to better correlations. This demonstrates that variability carries information through a given system. The concept of stochastic resonance suggests that there is an optimal level of noise for a given system after which system performance will degrade. We can see this conceptually in our model where the overall systembehavior is a balance between the noise introduced at each node and the starting variability. A system with more noise requires more initial variance to preserve information. In conclusion, we believe the intercellular variability is beneficial to the system, whereby it will improve the population response and decrease the sensitivity of the system to perturbation by other factors such as external stressors. Limitations of the study Sources of extrinsic variance include differences in gene/protein expression, unsynchronized cell cycles or circadian rhythms, differences in mitochondrial number, and environmental differences. Each of these po- tential sources (with the possible exception of mitochondrial content) result from or in differences in gene/ protein expression. For example, cell cycle phase is determined by the expression of proteins that regulate cell cycle kinases, and these ultimately trigger further changes in transcription and translation. Thus, it is difficult to identify a source of extrinsic variance that is ultimately not driven by differences in gene/protein expression that would meet the requisite criteria. These studies were performed by overexpressing key signaling components such as Akt. It is important to recognize that overexpression will have an effect on the system behavior. In this case we were not able to observe any direct effect of Akt expression levels on the behavior of GLUT4. Furthermore, the endogenous proteins are still present in these systems and their behavior has not been measured, as this is not feasible. Ideally the reporters would replace the endogenous protein. Finally, we acknowledge that in our system the expression level of the reporters adds an additional level of experimental noise and there will also be other sources of noise such as measurement noise. We have, where possible, accounted for these and do not believe these sources of noise change our conclusions. Resource availability Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, James Burchfield (james.burchfield@sydney.edu.au). ll OPEN ACCESS iScience 24, 102118, February 19, 2021 11 Materials availability Resources and reagents generated as a part of this study are available on request. Data and code availability All data and code are available on request. METHODS All methods can be found in the accompanying Transparent Methods supplemental file. SUPPLEMENTAL INFORMATION Supplemental Information can be found online at https://doi.org/10.1016/j.isci.2021.102118. ACKNOWLEDGMENTS This work was funded by Australian Research Council (ARC) project grant number DP180103482 awarded to D.E.J. and J.G.B., National Health and Medical Research Council (NHMRC) project grant number GNT1120201 awarded to D.E.J., and a Diabetes Australia research grant number Y19G-FAZD awarded to D.J.F. and J.G.B.. D.M.N. was supported by an Australian Postgraduate Award scholarship, and D.E.J. was supported by an NHMRC Senior Principal Research Fellowship. P.Y. and A.M.S. were supported by ARC Discovery Early Career Research Awards (DECRAs) (DE170100759, DE180101520), and P.Y. also by an NHMRC Investigator Grant (1173469). L.K.N. is supported by a Victorian Cancer Agency Mid-Career Research Fellowship (MCRF18026), an Investigator Initiated Research Scheme grant from National Breast Cancer Foundation (IIRS-20-094); and the Metcalf Venture Grants Scheme administered by Cancer Council Victoria, Australia. The authors acknowledge the facilities, and the scientific and technical assistance, of the Australian Microscopy and Microanalysis Research Facility at the Charles Perkins Centre, University of Syd- ney. The contents of the publishedmaterial are solely the responsibility of the individual authors and do not reflect the views of the funding agencies. AUTHOR CONTRIBUTION Conceptualization: D.N., J.G.B., D.E.J., and L.K.N; Methodology: D.N., P.Y., S.-Y.S., L.K.N., and J.G.B.; Formal analyses: D.N., P.Y., S.-Y.S., A.L.K., and J.G.B.; Investigation: D.N., P.Y., S.-Y.S., A.L.K., and J.G.B.; Data Curation: D.N., P.Y., S.-Y.S., and J.G.B.; Supervision: L.K.N., D.J., and J.G.B.; Writing – original draft: D.N. and J.G.B.; Writing – review and editing: D.N., P.Y., S.-Y.S., A.L.K., D.J.F., A.M.S., L.K.N., D.E.J., and J.G.B.; Supervision: L.K.N., J.G.B., and D.E.J.; Project administration: J.G.B.; Funding acquisition: D.E.J. DECLARATION OF INTERESTS The authors declare no competing interests. Received: November 16, 2020 Revised: January 7, 2021 Accepted: January 22, 2021 Published: February 19, 2021 REFERENCES Bai, L., Wang, Y., Fan, J., Chen, Y., Ji, W., Qu, A., Xu, P., James, D.E., and Xu, T. (2007). Dissecting multiple steps of GLUT4 trafficking and identifying the sites of insulin action. Cell Metab 5, 47–57. Balaban, N.Q., Merrin, J., Chait, R., Kowalik, L., and Leibler, S. (2004). Bacterial persistence as a phenotypic switch. Science 305, 1622–1625. Balbis, A., Baquiran, G., Bergeron, J.J.M., and Posner, B.I. (2000). Compartmentalization and insulin-induced translocations of insulin receptor substrates, phosphatidylinositol 3-kinase, and protein kinase B in rat liver. Endocrinology 141, 4041–4049. Biddy, B.A., Kong, W., Kamimoto, K., Guo, C., Waye, S.E., Sun, T., and Morris, S.A. (2018). Single-cell mapping of lineage and identity in direct reprogramming. Nature 564, 219–224. Brewer, P.D., Romenskaia, I., Kanow, M.A., and Mastick, C.C. (2011). Loss of AS160 Akt substrate causes Glut4 protein to accumulate in compartments that are primed for fusion in basal adipocytes. J. Biol. Chem. 286, 26287–26297. Buettner, F., Natarajan, K.N., Casale, F.P., Proserpio, V., Scialdone, A., Theis, F.J., Teichmann, S.A., Marioni, J.C., and Stegle, O. (2015). Computational analysis of cell-to-cell heterogeneity in single-cell RNA-sequencing data reveals hidden subpopulations of cells. Nat. Biotechnol. 33, 155–160. Burchfield, J.G., Lu, J., Fazakerley, D.J., Tan, S.-X., Ng, Y., Mele, K., Buckley, M.J., Han, W., Hughes, W.E., and James, D.E. (2013). Novel systems for dynamically assessing insulin action in live cells reveals heterogeneity in the insulin response. Traffic 14, 259–273. ll OPEN ACCESS 12 iScience 24, 102118, February 19, 2021 Cedersund, G., Roll, J., Ulfhielm, E., Danielsson, A., Tidefelt, H., and Stra˚lfors, P. (2008). Model- based hypothesis testing of key mechanisms in initial phase of insulin signaling. Plos Comput. Biol. 4, e1000096. Chelaru, M.I., and Dragoi, V. (2008). Efficient coding in heterogeneous neuronal populations. Proc. Natl. Acad. Sci. U S A 105, 16344–16349. Cheong, R., and Levchenko, A. (2010). Oscillatory signaling processes: the how, the why and the where. Curr. Opin. Genet. Dev. 20, 665–669. Chi, K.R. (2014). Singled out for sequencing. Nat. Methods 11, 13–17. Cline, G.W., Petersen, K.F., Krssak, M., Shen, J., Hundal, R.S., Trajanoski, Z., Inzucchi, S., Dresner, A., Rothman, D.L., and Shulman, G.I. (1999). Impaired glucose transport as a cause of decreased insulin-stimulated muscle glycogen synthesis in type 2 diabetes. N. Engl. J. Med. 341, 240–246. Cooper, S., and Bakal, C. (2017). Accelerating live single-cell signalling studies. Trends Biotechnol. 35, 422–433. Copps, K.D., and White, M.F. (2012). Regulation of insulin sensitivity by serine/threonine phosphorylation of insulin receptor substrate proteins IRS1 and IRS2. Diabetologia 55, 2565– 2582. Dibble, C.C., Asara, J.M., and Manning, B.D. (2009). Characterization of Rictor phosphorylation sites reveals direct regulation ofmTOR complex 2 by S6K1. Mol. Cell. Biol. 29, 5657–5670. Ebner, M., Lucic, I., Leonard, T.A., and Yudushkin, I. (2017). PI(3,4,5)P engagement restricts Akt activity to cellular membranes. Mol. Cell 65, 416– 431.e6. Evers, T.M.J., Hochane, M., Tans, S.J., Heeren, R.M.A., Semrau, S., Nemes, P., and Mashaghi, A. (2019). Deciphering metabolic heterogeneity by single-cell analysis. Anal. Chem. 91, 13314–13323. Fazakerley, D.J., Chaudhuri, R., Yang, P., Maghzal, G.J., Thomas, K.C., Krycer, J.R., Humphrey, S.J., Parker, B.L., Fisher-Wellman, K.H., Meoli, C.C., et al. (2018). Mitochondrial CoQ deficiency is a common driver of mitochondrial oxidants and insulin resistance. Elife 7, e23111. Gonzalez, E., Flier, E., Molle, D., Accili, D., and McGraw, T.E. (2011). Hyperinsulinemia leads to uncoupled insulin regulation of the GLUT4 glucose transporter and the FoxO1 transcription factor. Proc. Natl. Acad. Sci. U S A 108, 10162– 10167. Gonzalez, E., and McGraw, T.E. (2006). Insulin signaling diverges into Akt-dependent and -independent signals to regulate the recruitment/ docking and the fusion of GLUT4 vesicles to the plasma membrane. Mol. Biol. Cell 17, 4484–4493. Gorre, M.E., Mohammed, M., Ellwood, K., Hsu, N., Paquette, R., Rao, P.N., and Sawyers, C.L. (2001). Clinical resistance to STI-571 cancer therapy caused by BCR-ABL gene mutation or amplification. Science 293, 876–880. Greaves, M. (2015). Evolutionary determinants of cancer. Cancer Discov. 5, 806–820. Hoehn, K.L., Hohnen-Behrens, C., Cederberg, A., Wu, L.E., Turner, N., Yuasa, T., Ebina, Y., and James, D.E. (2008). IRS1-independent defects define major nodes of insulin resistance. Cell Metab. 7, 421–433. James, S.R., Downes, C.P., Gigg, R., Grove, S.J., Holmes, A.B., and Alessi, D.R. (1996). Specific binding of the Akt-1 protein kinase to phosphatidylinositol 3,4,5-trisphosphate without subsequent activation. Biochem. J. 315 (Pt 3), 709–713. Kashima, Y., Suzuki, A., Liu, Y., Hosokawa, M., Matsunaga, H., Shirai, M., Arikawa, K., Sugano, S., Kohno, T., Takeyama, H., et al. (2018). Combinatory use of distinct single-cell RNA-seq analytical platforms reveals the heterogeneous transcriptome response. Sci. Rep. 8, 3482. Keren, I., Shah, D., Spoering, A., Kaldalu, N., and Lewis, K. (2004). Specialized persister cells and the mechanism of multidrug tolerance in Escherichia coli. J. Bacteriol. 186, 8172–8180. Kim, C., Gao, R., Sei, E., Brandt, R., Hartman, J., Hatschek, T., Crosetto, N., Foukakis, T., and Navin, N.E. (2018). Chemoresistance evolution in triple-negative breast cancer delineated by single-cell sequencing. Cell 173, 879–893.e13. Klip, A., McGraw, T.E., and James, D.E. (2019). Thirty sweet years of GLUT4. J. Biol. Chem. 294, 11369–11381. Kohn, A.D., Summers, S.A., Birnbaum, M.J., and Roth, R.A. (1996). Expression of a constitutively active Akt Ser/Thr kinase in 3T3-L1 adipocytes stimulates glucose uptake and glucose transporter 4 translocation. J. Biol. Chem. 271, 31372–31378. Lasserre, R., Guo, X.-J., Conchonaud, F., Hamon, Y., Hawchar, O., Bernard, A.-M., Soudja, S.M., Lenne, P.-F., Rigneault, H., Olive, D., et al. (2008). Raft nanodomains contribute to Akt/PKB plasma membrane recruitment and activation. Nat. Chem. Biol. 4, 538–547. Lee, M.-C.W., Lopez-Diaz, F.J., Khan, S.Y., Tariq, M.A., Dayn, Y., Vaske, C.J., Radenbaugh, A.J., Kim, H.J., Emerson, B.M., and Pourmand, N. (2014). Single-cell analyses of transcriptional heterogeneity during drug tolerance transition in cancer cells by RNA sequencing. Proc. Natl. Acad. Sci. 111, E4726–E4735. Lizunov, V.A., Stenkula, K.G., Blank, P.S., Troy, A., Lee, J.-P., Skarulis, M.C., Cushman, S.W., and Zimmerberg, J. (2015). Human adipose cells in vitro are either refractory or responsive to insulin, reflecting host metabolic state. PLoS One 10, e0119291. McDonnell, M.D., and Abbott, D. (2009). What is stochastic resonance? Definitions, misconceptions, debates, and its relevance to biology. PLoS Comput. Biol. 5, e1000348. Miragaia, R.J., Gomes, T., Chomka, A., Jardine, L., Riedel, A., Hegazy, A.N.,Whibley, N., Tucci, A., Chen, X., Lindeman, I., et al. (2019). Single-cell transcriptomics of regulatory T cells reveals trajectories of tissue adaptation. Immunity 50, 493–504.e7. Navin, N., Kendall, J., Troge, J., Andrews, P., Rodgers, L., McIndoo, J., Cook, K., Stepansky, A., Levy, D., Esposito, D., et al. (2011). Tumour evolution inferred by single-cell sequencing. Nature 472, 90–94. Ng, Y., Ramm, G., Lopez, J.A., and James, D.E. (2008). Rapid activation of Akt2 is sufficient to stimulate GLUT4 translocation in 3T3-L1 adipocytes. Cell Metab. 7, 348–356. Nolan-Stevaux, O., Tedesco, D., Ragan, S., Makhanov, M., Chenchik, A., Ruefli-Brasse, A., Quon, K., and Kassner, P.D. (2013). Measurement of cancer cell growth heterogeneity through lentiviral barcoding identifies clonal dominance as a characteristic of in vivo tumor engraftment. PLoS One 8, e67316. Norris, D.M., Yang, P., Krycer, J.R., Fazakerley, D.J., James, D.E., and Burchfield, J.G. (2017b). An improved Akt reporter reveals intra- and inter- cellular heterogeneity and oscillations in signal transduction. J. Cell Sci. 130, 2757–2766. Park, W.S., Heo, W.D., Whalen, J.H., O’Rourke, N.A., Bryan, H.M., Meyer, T., and Teruel, M.N. (2008). Comprehensive identification of PIP3- regulated PH domains from C. elegans to H. sapiens by model prediction and live imaging. Mol. Cell 30, 381–392. Petitpre´, C., Wu, H., Sharma, A., Tokarska, A., Fontanet, P., Wang, Y., Helmbacher, F., Yackle, K., Silberberg, G., Hadjab, S., and Lallemend, F. (2018). Neuronal heterogeneity and stereotyped connectivity in the auditory afferent system. Nat. Commun. 9, 3691. Shalek, A.K., Satija, R., Adiconis, X., Gertner, R.S., Gaublomme, J.T., Raychowdhury, R., Schwartz, S., Yosef, N., Malboeuf, C., Lu, D., et al. (2013). Single-cell transcriptomics reveals bimodality in expression and splicing in immune cells. Nature 498, 236–240. Sharick, J.T., Jeffery, J.J., Karim, M.R., Walsh, C.M., Esbona, K., Cook, R.S., and Skala, M.C. (2019). Cellular metabolic heterogeneity in vivo is recapitulated in tumor organoids. Neoplasia 21, 615–626. Sto¨ckli, J., Davey, J.R., Hohnen-Behrens, C., Xu, A., James, D.E., and Ramm, G. (2008). Regulation of glucose transporter 4 translocation by the Rab guanosine triphosphatase-activating protein AS160/TBC1D4: role of phosphorylation and membrane association. Mol. Endocrinol. 22, 2703–2715. Su, Z., Burchfield, J.G., Yang, P., Humphrey, S.J., Yang, G., Francis, D., Yasmin, S., Shin, S.-Y., Norris, D.M., Kearney, A.L., et al. (2019). Global redox proteome and phosphoproteome analysis reveals redox switch in Akt. Nat. Commun. 10, 5486. Tabula Muris Consortium (2018). Overall coordination, Logistical coordination, Organ collection and processing, Library preparation and sequencing, Computational data analysis, Cell type annotation, Writing group, Supplemental text writing group, Principal investigators,. Nature 562, 367–372, Single-cell transcriptomics of 20 mouse organs creates a Tabula Muris. Tan, S.-X., Fisher-Wellman, K.H., Fazakerley, D.J., Ng, Y., Pant, H., Li, J., Meoli, C.C., Coster, A.C.F., Sto¨ckli, J., and James, D.E. (2015). Selective insulin resistance in adipocytes. J. Biol. Chem. 290, 11337–11348. ll OPEN ACCESS iScience 24, 102118, February 19, 2021 13 Tan, S.-X., Ng, Y., Burchfield, J.G., Ramm, G., Lambright, D.G., Sto¨ckli, J., and James, D.E. (2012a). The Rab GTPase-activating protein TBC1D4/AS160 contains an atypical phosphotyrosine-binding domain that interacts with plasma membrane phospholipids to facilitate GLUT4 trafficking in adipocytes. Mol. Cell. Biol. 32, 4946–4959. Tan, S.-X., Ng, Y., Meoli, C.C., Kumar, A., Khoo, P.-S., Fazakerley, D.J., Junutula, J.R., Vali, S., James, D.E., and Sto¨ckli, J. (2012b). Amplification and demultiplexing in insulin-regulated Akt protein kinase pathway in adipocytes. J. Biol. Chem. 287, 6128–6138. Tremblay, F., Lavigne, C., Jacques, H., and Marette, A. (2001). Defective insulin-induced GLUT4 translocation in skeletal muscle of high fat-fed rats is associated with alterations in both Akt/protein kinase B and atypical protein kinase C (zeta/lambda) activities. Diabetes 50, 1901– 1910. Wada, T., Hironaka, K.-I., Wataya, M., Fujii, M., Eto, M., Uda, S., Hoshino, D., Kunida, K., Inoue, H., Kubota, H., et al. (2020). Single-cell information analysis reveals that skeletal muscles incorporate cell-to-cell variability as information not noise. Cell Rep. 32, 108051. Wang, D., and Bodovitz, S. (2010). Single cell analysis: the new frontier in ‘‘omics. Trends Biotechnol. 28, 281–290. Wang, Y., Waters, J., Leung, M.L., Unruh, A., Roh, W., Shi, X., Chen, K., Scheet, P., Vattathil, S., Liang, H., et al. (2014). Clonal evolution in breast cancer revealed by single nucleus genome sequencing. Nature 512, 155–160. Xiong, W., Jordens, I., Gonzalez, E., and McGraw, T.E. (2010). GLUT4 is sorted to vesicles whose accumulation beneath and insertion into the plasma membrane are differentially regulated by insulin and selectively affected by insulin resistance. Mol. Biol. Cell 21, 1375–1386. Yang, P., Humphrey, S.J., Cinghu, S., Pathania, R., Oldfield, A.J., Kumar, D., Perera, D., Yang, J.Y.H., James, D.E., Mann, M., and Jothi, R. (2019). Multi- omic profiling reveals dynamics of the phased progression of pluripotency. Cell Syst. 8, 427– 445.e10. Zhang, Z., Milias-Argeitis, A., and Heinemann, M. (2018). Dynamic single-cell NAD(P)H measurement reveals oscillatory metabolism throughout the E. coli cell division cycle. Sci. Rep. 8, 2162. Zhao, J., Zong, W., Zhao, Y., Gou, D., Liang, S., Shen, J., Wu, Y., Zheng, X., Wu, R., Wang, X., et al. (2019). In vivo imaging of b-cell function reveals glucose-mediated heterogeneity of b-cell functional development. Elife 8, e41540. ll OPEN ACCESS 14 iScience 24, 102118, February 19, 2021 iScience, Volume 24 Supplemental Information Signaling Heterogeneity is Defined by Pathway Architecture and Intercellular Variability in Protein Expression Dougall Norris, Pengyi Yang, Sung-Young Shin, Alison L. Kearney, Hani Jieun Kim, Thomas Geddes, Alistair M. Senior, Daniel J. Fazakerley, Lan K. Nguyen, David E. James, and James G. Burchfield Supplemental Information Supplemental Figures and Legends Supplementary Figure S1. Single cell gene expression data from 3T3-L1 fibroblasts reveals that the expression of individual genes can vary by more than an order of magnitude between cells. Related to Figure 5. A)Histograms showing the magnitude of expression differences of genes between cells. Expression of each gene was log2 transformed and the magnitude of difference for each gene was calculated between the minimum and maximum value, the top and bottom 10th percentile and the top and bottom quartile.B) The min:max fold change verse the average gene expression for a subset of genes associated with insulin signalling. Supplementary Figure S2. Example response profiles and correlations from the generalised multi-tier signalling networks. Related to Figure 6A and B. A) and B) are without feedback, C) and D) with feedback. Supplementary Figure S3. Negative feedback induces an overshoot, suppresses the overal signal and decreases the response time. The strength of feedback in a representative dataset with fixed expression levels as varied from weak to strong. Related to Figure 6. The response at each node are shown normalised to the absolute max within each node (upper panels) or to the max of the individual curve (lower panels). Supplementary Figure S4. The relative variability (coefficient of variation; CV) in the steady state response to 1nM insulin is greater for GLUT4 than Akt. Related to Figures 3 and 6. Points represent the CV for each each experiment presented in Figure 2. Bar is the mean +/- the SD. *** p<0.001 by unpaired t-test with Welch’s correction. Supplementary Figure S5. Increased variance in the primary node (S1) leads to better information transmission through signalling cascades. Related to Figure 7. A) Scatter plots showing correlation between Rho (correlation between aS1 and downstream nodes) and coefficient of variance for aS1, for a single parameter set where the expression variation of the primary node was gradually increased up to 1000 folds n the negative feedback model. The expression of network nodes and phosphatases were allowed to randomly vary within a maximum of 5 folds of nominal values (n = 50). Note that 20 folds of variation of the primary node generated very small CV of aS1.B)Alternative representation of correlation and CV for gradual increase of the primary node variation as in (A) for non-feedback system. Transparent Methods Cell culture 3T3-L1 murine fibroblasts sourced from the American Type Culture Collection (ATCC, Manassas, VA) were maintained as described previously (Norris et al., 2017a) and routinely tested and proven to be free of contamination. In brief, cells were cultured in high-glucose Dulbecco’s modified Eagle’s medium (DMEM), supplemented with 10% fetal bovine serum (FBS) and Glutamax (sourced from Gibco) at 37°C in the presence of 10% CO2. Differentiation of confluent fibroblasts into adipocytes was induced with the addition of 350 nM insulin, 250 nM dexamethasone, 0.5 mM 3-isobutyl-1-methylxanthine and 400 nM biotin for 3 days (all of which were from Sigma-Aldrich). Cells were then incubated in media containing 350 nM insulin for a further 3 d, and refreshed with 10% FBS, DMEM and glutamax every 2 d afterwards. Experiments were performed from day 8-10 post-differentiation. Immunofluorescence Cells were seeded into 8 well Ibidi u-slides, coated with matrigel, seven d post-differentiation and then allowed to attach for 2 days. 9 days post-differentiation, cells were serum-starved in DMEM (with 0.2% BSA and Glutamax) for 2 hours at 37°C with 10% CO2 and then stimulated with 100 nM insulin. The coverslips were then briefly immersed in an ice-cold PBS bath and instantly fixed with 4% paraformaldehyde at room temperature for 15 min. Cells were then washed twice with room temperature PBS and quenched with 100 mM Glycine for 10 min. The cells were washed twice more with room temperature PBS, and incubated in blocking and permeabilising buffer (PBS with 2% BSA and 0.1% saponin) for 30 min. Cells were then incubated overnight at 4°C in blocking and permeabilising buffer containing the anti-phospho-T308-Akt primary antibody (rabbit; sourced from Cell Signalling Technologies - #9275) at a dilution of 1:100. The following day, the cells were washed with blocking and permeabilising buffer 5 times, and then incubated with anti-rabbit Alexa-488 (1:500) at room temperature for 30 min in the dark. Cells were washed 5 more times with PBS and then stored and imaged in PBS, 5% glycerol and 2.5% Dabco. Live cell imaging For live-cell imaging experiments, plasmids were delivered by electroporation 7 d post-differentiation(see (Norris et al., 2017b)). Day 7 3T3-L1 adipocytes were trypsinized, pelleted at 150 x g for 5 min, and washed three times with PBS before resuspension in electroporation solution (20 mMHEPES, 135 mMKCl, 2 mMMgCl2, 0.5%(w/v) Ficoll 400, 1% (v/v) DMSO, 2 mMATP, and 5 mMGSH,pH 7.6) and 5µg of plasmid DNA. Cells were then electroporated at 200 mV for 20 ms and plated onto Matri-gel-coated 35-mm Ibidi glass-bottom m-Dishes and seeded onto Matrigel coated 35 mm μ-dishes (Ibidi). After24 h, cells were incubated with basal medium (DMEM without FBS) for 2 h. Fol-lowing this media was replaced with KRP+buffer (KRP, 10 mM glucose andessential amino acids) and dishes were placed onto the stage of a Nikon TiEmicroscope equipped with an OKOlab microscope enclosure maintained at 37 °C.TIRF was achieved using a Nikon hTIRF module. TagRFP-T was stimulated with a568 nm laser angled at 71 °C and emission was captured on an Andor 888 emCCDcamera after passing through a 610/50 nmfilter. Buffer switching was performedusing a customfluidic setup. and seeded on matrigel coated coverslips. Cells were serum starved for 1.5 - 2 h in DMEM supplemented with 0.2% BSA and Glutamax and then transferred to a modified Krebs Ringer phosphate buffer (KRP) buffer consisting of 120 mM NaCl, 0.6 mM Na2HPO4, 0.4 mM NaH2PO4, 6 mM NaCl, 1.2 mM MgSO4, 12.5 mM HEPES, 1 mM CaCl2, 10 mM glucose, 0.2 % w/v BSA and Glutamax at pH 7.4 for imaging. Assessment of live signalling responses were performed at 37°C in the absence of CO2. Imaging and computational analysis Images were analysed using FIJI (Schindelin et al., 2012). All downstream computational analyses were conducted using R programming environment. First, quantifications of GLUT4 and Akt across the time- crouse were converted to ratios by dividing the quantification at each time point with the mean of those that prior to the treatment of insulin. Next, to combine GLUT4 and Akt, time-course profiles generated from multiple biological independent experiments, a locally weighted scatterplot smoothing (LOESS) model (Cleveland and Devlin, 1988) was fitted to each individual cell to interpolate and align time points across all experiments. Principal component analysis (PCA) was subsequently conducted to identify if there were batch efforts between individual experiments. Specifically, we summarised response profiles of GLUT4 and Akt in cells across measured time points into principal components (PCs) and investigated if the cells are grouped by their experimental batches in the PC space. For correlation analysis, we extracted the maximum response (max), time taken for activation (lag, defined as time taken to reach 10% of max), and half time (defined as time taken to reach 50% of max) from the GLUT4 profile of each cell. Similarly, we extracted the max, the area under the curve of response profiles (AUC) and the shape of the response profiles (min-max normalised AUC) from the Akt profile of each cell. Fuzzy c-means clustering was used to analyse the shape of the Akt response profiles and two clusters with distinctive temporal shapes were discovered. The importance of Akt characteristics are also assessed using a random forest regression model where they are ranked by their predictiveness to GLUT4 max. For model inference, four models were defined with GLUT4 max as response variable and various combinations of Akt profile features as explanatory variables. These models were subsequently fitted to the data to identify explanatory variables and their combinations that were significantly associated with GLUT4 max response. Single cell data and calculation of fold change The processed single-cell expression matrix was downloaded from Tabula Muris Consortium (https://doi.org/10.1038/s41586-018-0590-4; https://tabula-muris.ds.czbiohub.org/). Specifically, gene counts derived from SMART-Seq2 RNAseq technology were used. The expression matrix was normalised using `logNormCounts` function from the `scater` package (McCarthy et al., 2017). We used cell type annotations defined from the original paper (Tabula Muris Consortium et al., 2018). For the current analysis, cell types with more than 600 cells were used, which resulted in a total of 18 cell types resolved by tissue origin. For each cell type, we calculated the intercellular log2 fold change in expression by calculating the difference in the median of gene expression from the top and bottom 25% quartiles of non-zero expression values. The log2 fold change in gene expression was calculated for each gene. Genes annotated to be part of the GO:0007169 pathway were extracted using `org.Mm.eg.db` (Marc Carlson (2020). org.Mm.eg.db: Genome wide annotation for Mouse. R package version 3.11.4.) and `GO.db` packages (Marc Carlson (2020). GO.db: A set of annotation maps describing the entire Gene Ontology. R package version 3.11.4.). For the genes in the gene set, the cell type specific intracellular log2 fold change expression values were averaged to generate an average log2 fold change across all cell types. The average gene expression level and coefficient of variation was calculated in the same way. Kinetic modelling of generalized multi-tier signalling networks To computationally analyse the dynamics of a generalized signalling cascade, we constructed a kinetic model of a multi-tier, phosphorylation cascade using ordinary differential equations (ODEs) that are formulated based on Michaelis-Menten kinetics. The model’s schematic diagram containing the model reactions is given in Figure 5A and C. Model description including ODEs and rate equations used for simulations are given in Supplementary Tables S1-4. The model was implemented and numerically simulated in MATLAB ® (The MathWorks. Inc. 2020a) and IQM Tools (http://www.intiquan.com/intiquan- tools/). IQM Tools was used to convert a standard IQM model into C code, which is compiled to an executable MEX file for model simulation. For the simulation analysis we generated random kinetic parameter sets (n=50) within wide ranges: kc = 0.001 – 100, Vm = 0.1 – 1000, and Km = 0.1 – 1000, and the expression of network nodes (Si = 100, i=1...6) and Phosphatase (Pi = 100, i=1...6) was allowed to randomly vary within a maximum of 5 folds of nominal values. We assumed that the cascade is stimulated by an input signal (in this case 10nM insulin) in all time-course simulations. Supplemental References McCarthy, D.J., Campbell, K.R., Lun, A.T.L., Wills, Q.F., 2017. Scater: pre-processing, quality control, normalization and visualization of single-cell RNA-seq data in R. Bioinformatics 33, 1179–1186. Norris, D.M., Geddes, T.A., James, D.E., Fazakerley, D.J., Burchfield, J.G., 2017a. Glucose Transport: Methods for Interrogating GLUT4 Trafficking in Adipocytes, in: Methods in Molecular Biology. pp. 193– 215. Norris, D.M., Geddes, T.A., James, D.E., Fazakerley, D.J., Burchfield, J.G., 2017b. Glucose Transport: Methods for Interrogating GLUT4 Trafficking in Adipocytes, in: Methods in Molecular Biology. pp. 193– 215. Schindelin, J., Arganda-Carreras, I., Frise, E., Kaynig, V., Longair, M., Pietzsch, T., Preibisch, S., Rueden, C., Saalfeld, S., Schmid, B., Tinevez, J.-Y., White, D.J., Hartenstein, V., Eliceiri, K., Tomancak, P., Cardona, A., 2012. Fiji: an open-source platform for biological-image analysis. Nat. Methods 9, 676–682. Tabula Muris Consortium, Overall coordination, Logistical coordination, Organ collection and processing, Library preparation and sequencing, Computational data analysis, Cell type annotation, Writing group, Supplemental text writing group, Principal investigators, 2018. Single-cell transcriptomics of 20 mouse organs creates a Tabula Muris. Nature 562, 367–372.