Influence of node abundance on signaling network state and dynamics analyzed by mass cytometry.

Signaling networks are key regulators of cellular function. Although the concentrations of signaling proteins are perturbed in disease states, such as cancer, and are modulated by drug therapies, our understanding of how such changes shape the properties of signaling networks is limited. Here we couple mass-cytometry-based single-cell analysis with overexpression of tagged signaling proteins to study the dependence of signaling relationships and dynamics on protein node abundance. Focusing on the epidermal growth factor receptor (EGFR) signaling network in HEK293T cells, we analyze 20 signaling proteins during a 1-h EGF stimulation time course using a panel of 35 antibodies. Data analysis with BP-R2, a measure that quantifies complex signaling relationships, reveals abundance-dependent network states and identifies novel signaling relationships. Further, we show that upstream signaling proteins have abundance-dependent effects on downstream signaling dynamics. Our approach elucidates the influence of node abundance on signal transduction networks and will further our understanding of signaling in health and disease.

Signaling networks are at the core of cellular information processing and transform external signals into cellular responses. Signals are transduced by modulating enzymatic activities mainly via protein phosphorylation, and cells implement sophisticated mechanisms, such as feedback loops, pathway crosstalk, and differential enzyme localization, to integrate signals and drive cellular processes and physiological outputs. The abundance of individual signaling pathway components (nodes) is central to the activity and output of a signaling network1. Changes in node abundance are tightly regulated and control biological programs such as stem cell differentiation and embryogenesis2. Abundance deregulation of particular signaling network nodes via genomic, transcriptional, or post-transcriptional regulatory defects3–5 underlies human diseases, the prime example being cancer6. Copy number alterations of genes encoding critical proteins7–9, independent of mutations that constitutively change enzymatic activity10, drive progression of many cancer types. Genomic instability in cancer cells causes abnormally broad distributions of signaling protein abundances in a given tumor11, yet the consequences of the protein abundance levels on signaling properties is poorly understood limiting our ability to rationally design therapies.

The epidermal growth factor receptor (EGFR) signaling network is affected by gene copy number alterations that deregulate protein abundances (e.g., of EGFR, HER2, ERK and AKT) in a number of cancer types7–9. EGFR signaling controls cell growth, motility, survival, differentiation, and metabolism12. Many drugs target the activity of the EGFR signaling network13,14. The receptor tyrosine kinase (RTK) function of EGFR is activated by its dimerization upon ligand binding. EGFR auto-phosphorylation recruits adaptor proteins that typically activate the MAPK/ERK and AKT signaling pathways. The MAPK/ERK branch activates the GTPase RAS, which triggers a kinase phosphorylation cascade consisting of RAF, MEK, ERK, and p90RSK. The output of the MAPK/ERK branch is transcription of genes regulating growth and division15,16. Signal transduction through the AKT branch starts by PI3K activation, producing PIP3, which recruits AKT and PDK1 to the plasma membrane. PDK1 phosphorylates AKT15,17, which mediates signaling through the mTORC1 complex to modulate translation via p70S6K and 4EBP117. Other AKT targets are GSK3β, PRAS40, and TSC2. The AKT pathway controls cell survival, proliferation, and migration17. STAT proteins and the PKC pathway can also be activated by EGFR-mediated signaling18,19. EGFR signaling involves crosstalk and feedback loops both internally (e.g., active ERK attenuates upstream RAF or MEK signaling via negative feedback)15 and with other signaling pathways (e.g., WNT and TGF-β pathways)20,21.

Classically, two approaches are used to characterize the effect of proteins on signal transduction. The first approach analyzes cell populations. Here, western blotting, mass spectrometry, RNA-microarrays, and synthetic lethality screens are used to identify signaling relationships22–24. Protein-protein interaction analyses are used to determine which proteins in a network directly interact23,25. Population-based methods yield a comprehensive view of signaling but are difficult to use in analysis of protein abundance dependencies due to inherent limitations: Proteins must be expressed at different abundances or cells must be sorted to yield a non-continuous abundance titration. Such methods result in a large number of samples and cell-to-cell protein abundance variations within each sample remain masked. The second approach studies signaling relationships in single cells. Here fluorescence microscopy and flow cytometry (FACS) are used with a variety of assays, including proximity ligation assay (PLA)26 or fluorescence resonance energy transfer (FRET)27. These approaches allow study of signaling relationships and dynamics through time and space; however, only a few signaling nodes can be measured simultaneously.

A recently developed single-cell analysis technology, called mass cytometry, allows for the simultaneous measurement of over 40 signaling nodes in single cells using metal-isotope tagged antibodies28,29. This capability makes mass cytometry uniquely suited to comprehensively query the function of nodes in signaling networks within heterogeneous cell populations. Mass cytometry is quantitative and, in combination with mass-tag cellular barcoding (MCB), a powerful screening tool28. Algorithms to analyze multiplexed single-cell mass cytometry data allow quantification of signaling relationships, therefore helping to decipher the highly complex network behaviors that operate even in simple biological systems30.

Here, we coupled protein overexpression with mass cytometry to measure the effect of varying node abundance on the activation state and signaling relationships of an unstimulated EGFR signaling network, as well as the signaling dynamics of the network in response to EGF stimulation. We exploited the finding that transient protein overexpression in a cell population typically produces a continuous abundance range of the target protein over four orders of magnitude. We overexpressed 20 central EGFR signaling network proteins individually in human embryonic kidney (HEK) 293T cells, sampled during an EGF stimulation time course over 60 minutes totaling 360 conditions. An average of 11,000 cells per condition was analyzed with a panel of 35 antibodies to provide a comprehensive single-cell proteomic EGFR network analysis. To identify signaling relationships in this dataset, we developed a statistical measure that we call 'binned pseudo R-squared' (BP-R2) that recapitulated known signaling relationships and identified relationships that were –to the best of our knowledge- not described previously. Thus, our experimental and computational approach enables study of how the strength and dynamics of signal transduction are tuned by node abundances.

Results

Analyzing continuous protein abundance dependencies

To systematically identify and characterize protein abundance-dependent signaling relationships, dynamics, and network activation states, we exploited the variation and large dynamic range of protein abundance induced by transient transfection and used mass cytometry to quantify the abundance of the transfected protein of interest (POI) in conjunction with comprehensive signaling network readouts in single cells. We cloned POIs genes into vectors containing a cytomegalovirus (CMV) promoter and a GFP-tag sequence31 to transiently overexpress GFP-tagged POIs in HEK293T cells (Fig. 1a). The tagged protein abundance was measured by mass cytometry using an anti-GFP antibody (Fig. 1a). Ordering the measured cells based on the GFP signal provided a continuous POI titration (Fig. 1b). Typically, not all cells were transfected, yielding an internal control for every experiment. To measure the single-cell EGFR signaling network states, we designed and validated a panel of 35 antibodies that mostly detect phosphorylation sites on signaling proteins (Supplementary Tables 1-3). These data were used to determine the abundance dependencies of network activation state and signaling dynamics (Fig. 1b).

Figure 1

Workflow of abundance-dependent network analysis. (a) Experimental workflow. Signaling POIs are cloned into vectors containing a CMV promoter and a GFP-tag sequence to transiently overexpress GFP-tagged POIs in HEK293T cells. We quantify anti-GFP antibody as readout of POI-GFP abundance, together with other 35 markers, by mass cytometry. (b) Data analysis workflow. Cells were ordered based on the GFP signal, providing a continuous POI titration, which was then coupled to other signaling markers to determine the abundance dependencies of network activation state and signaling dynamics in the network after transfection. The network in the illustration does not represent an actual biological example.

To validate our system we confirmed that, first, the GFP tag was reliably detected by mass cytometry (Supplementary Fig. 1); second, the GFP tag did not affect the localization and activity of the POI (Supplementary Fig. 2, 3, Supplementary Table 4, Supplementary File 1); third, POI expression levels were linearly related to GFP abundance, validating GFP as readout of the total POI abundance (Supplementary Fig. 4a, c); fourth, POI overexpression for 18 hours (i.e., the time point of our experiments) did not alter the underlying network structure (Supplementary Fig. 4b, c); fifth, the antibody-based GFP quantification by mass cytometry was comparable to FACS (Supplementary Fig. 5); sixth, the cell culture media and cell detachment did not alter signaling processing in the EGFR network (Supplementary Fig. 6, 7); and, seventh, the levels of the GFP-tagged POIs were stable during the 1-hour EGF stimulation time course (Supplementary Fig. 8, Supplementary Video 1). We also found that the method is robust and highly reproducible as evidenced by the high concordance between the three individual experiment replicates (Supplementary Fig. 9, Supplementary File 2).

KRASG12V and MEK1DD abundance effect on signaling

We first studied a well-known signaling circuit: Constitutively active mutants of KRAS and MEK1 (KRASG12V and MEK1DD) lead to ERK phosphorylation and activate components downstream in the MAPK/ERK pathway. As expected, we found that overexpression of KRASG12V-GFP or MEK1DD-GFP increased phosphorylation on Thr202 and Tyr204 of ERK1/2 (Fig. 2a). Our approach also elucidated the abundance-dependent effects on these signaling relationships: The relationship between KRASG12V-GFP and p-ERK1/2 was bow-like as high levels of KRASG12V-GFP corresponded to reduced phosphorylation of ERK1/2. By contrast, the MEK1DD-GFP abundance relationship with p-ERK1/2 was monotonic as p-ERK1/2 increased with MEK1DD-GFP expression (Fig. 2a). These results verified the oncogenic activation of p-ERK1/2 induced by KRASG12V and MEK1DD.

Figure 2

MAPK/ERK pathway mutants induce oncogenic signaling. (a) Biaxial plots of GFP, representing the abundance of the overexpressed mutant POIs, versus abundance of phosphorylation on Thr202/Tyr204 on ERK1/2. Constitutively active KRASG12V-GFP shows a downregulation on Thr202/Tyr204 on ERK1/2 at the highest levels of KRASG12V-GFP. Constitutively active MEK1DD-GFP directly phosphorylates Thr202/Tyr204 on ERK1/2, and the abundance of the POI-GFP is correlated with amount of ERK1/2 phosphorylated at these sites. The FLAG-GFP control does not affect ERK phosphorylation sites. (b) The abundances of measured phosphorylation sites are plotted over the range of the KRASG12V-GFP and MEK1DD-GFP expression. Phosphorylation sites of the same pathway (e.g., on ERK1/2 and p90RSK, AKT and GSK3β, or p38 and JNK) show similar trends. An individual experiment is shown here. Plots for 3 replicates are shown in Supplementary Fig. 9b-e. (c) Strong single-cell correlations within biaxial plots indicate co-regulated phosphorylation sites. (d) Unchanged and reduced correlations indicate unrelated phosphorylation sites. (e) and (f) Heat maps showing for all pairs of measured markers the change in Fisher-transformed Spearman correlation values for overexpression of (e) KRASG12V-GFP and (f) MEK1DD-GFP when compared to the FLAG-GFP overexpression control. (g) and (h) BP-R2 scores and Spearman correlations of bin medians for all measured markers in cells where (g) KRASG12V-GFP or (h) MEK1DD-GFP was overexpressed overlaid on a literature-based graph of canonical signaling pathways14,15,21,23,44,45,35,46–48. Strong relationships identified from the BP-R2 analysis are plotted on the signaling maps as colored circles. The sizes of circles indicate relationship strengths quantified by BP-R2. The directionalities of relationships, as judged by Spearman correlation of bin medians, are shown by the color of the circles (positive correlation indicates that cells show generally increasing marker levels, and a negative correlation indicates decreasing marker levels as POI-GFP levels increase). For (e) to (h), data from 3 individual experiment replicates were used.

Next, we analyzed the impact of KRASG12V-GFP and MEK1DD-GFP abundance on all measured phosphorylation sites. We divided the measured cells into 10 bins according to the GFP signals and plotted the bin medians (Fig. 2b, Supplementary Fig. 9b-e). This analysis revealed that the phosphorylation site abundances on ERK1/2 and its direct downstream target Ser380 of p90RSK had similar relationships to the abundances of KRASG12V-GFP or MEK1DD-GFP. Phosphorylation of AKT on Ser473 and its direct target Ser9 of GSK3β also had parallel trends and showed reduced levels when the MAPK/ERK signal peaked, suggesting inter-pathway regulation. We also observed increased JNK phosphorylation on Thr183/Tyr185 induced by the KRASG12V mutant (Fig. 2b) as reported previously32. This shows that our approach recapitulates known signaling relationships and identifies abundance-determined signaling responses.

We then systematically evaluated signaling relationships between all pairs of measured markers modulated by KRASG12V-GFP or MEK1DD-GFP overexpression. We exploited the fact that overexpression of one protein increases signaling (i.e., phosphorylation levels) and thus expands the dynamic range of many measured markers (Fig. 2c). This enabled the use of correlation analysis to distinguish signaling relationships (high correlation) from biological and technical noise (low correlation). For example, overexpression of KRASG12V-GFP resulted in an increased Spearman correlation between p-ERK1/2 and p-p90RSK compared to control (Fig. 2c), whereas ERK-independent phosphorylation sites, such as Tyr551 of BTK/ITK, showed low correlation with p-ERK1/2 levels in both control and overexpression conditions (Fig. 2d).

Identifying changes in pairwise Spearman correlations for all measured markers in the KRASG12V-GFP and MEK1DD-GFP overexpression data compared to the FLAG-GFP control enabled systematic analysis of signaling relationship patterns (Fig. 2e, f). Phosphorylation levels of proteins in the MAPK/ERK pathways showed strong increases in correlation, and pathway members clustered together (Fig. 2e, f, green squares). We also observed that phosphorylations of MAPK/p38 pathway members and STAT proteins (STAT1 and STAT5) were increasingly correlated with levels of MAPK/ERK pathway members as MEKDD-GFP levels increased (Fig. 2f, purple rectangle), indicating crosstalk between MAPK and STAT pathways. These results reveal relationships among many measured markers and show that increases in correlation reflect pathways and grouped biological processes.

Automated analysis of abundance-induced signaling

Spearman correlation analysis can uncover strictly monotonic relationships between phosphorylation levels on signaling proteins; however, protein abundance-dependent signaling responses can be complex (Fig. 2a, see KRASG12V). We therefore developed a density-independent measure termed 'binned pseudo R-squared' (BP-R2) to quantify the strengths of relationships between the abundance of a POI and measured phosphorylation sites. BP-R2 creates 10 bins across the POI-GFP expression range and calculates the relationship strength considering bin medians and the global mean (Supplementary Fig. 10a, b, Methods, Supplementary Software). Using the BP-R2 values for all negative controls, a cutoff for strong signaling relationships was determined (Supplementary Fig. 10c). Benchmarking BP-R2 in identifying strong signaling relationships from the overexpression datasets showed that BP-R2 outperformed methods often used for this task30,33 (Supplementary Fig. 11a, b). The strong relationships identified by BP-R2 were plotted in a two-dimensional layout guided by canonical pathways (Fig. 2g, h). The directionality of measured signaling relationships was determined by Spearman correlation of the bin medians (Supplementary Fig. 10b, Methods). A positive correlation indicates that cells show generally increasing marker levels and a negative correlation indicates generally decreasing marker levels as POI-GFP levels increase.

Analysis of KRASG12V-GFP and MEK1DD-GFP overexpression versus all measured markers using BP-R2 revealed strong, positively correlated relationships of MEKDD-GFP to downstream MAPK/ERK pathway nodes. KRASG12V-GFP levels, although also positively correlated with MAPK/ERK nodes, exhibited the same, but weaker relationships (Fig. 2a, b, g, h). Together, these results suggest that feedback regulation of upstream MAPK nodes differs between the studied mutants. Additionally, this network view revealed that MEK1DD-GFP abundance had a strong positive impact on nodes in the MAPK/p38 pathway; the previously observed KRASG12V-induced phosphorylation of JNK32 was dependent on KRASG12V abundance (Fig. 2g, h). These results show that overexpression of signaling proteins, in conjunction with BP-R2 and correlation analysis, identifies known relationships and is a valid platform for discovery of signaling relationships in a comprehensive and abundance-dependent manner.

Node abundance dependency analyses of the EGFR network

To study the node abundance dependency of signaling relationships and dynamics in the EGFR signaling network, we overexpressed 20 EGFR-related signaling proteins individually in HEK293T cells (Table 1). Each of the 20 GFP-tagged POIs was validated in previous studies (Supplementary Table 5) and in our system (Supplementary Fig. 2, 3, Supplementary File 1). 18 hours after transfection, we treated cells with EGF and quantified signaling by mass cytometry over a 60-min time course. To exclude signaling relationships caused by channel-to-channel spillover, we applied a stringent experimental filter (Supplementary Fig. 12, Methods). The median marker intensities during the time course are shown in Supplementary Fig. 13a. Based on these data we performed two sets of analyses. In the first, we used BP-R2 analysis and Spearman correlations to evaluate how the abundance of overexpressed proteins influenced phosphorylation at the measured sites (Fig. 3, Supplementary Fig. 13b and Supplementary Files 2-4). In the second, we examined how features of signaling dynamics depend on protein abundance (Fig. 4).

Table 1

Overexpressed signaling proteins

Overexpressed proteins	Gene ID	UniProt Entry
SRC	SRC	P12931
PDK1	PDPK1	O15530
AKT1	AKT1	P31749
GSK3β	GSK3B	P49841
MKK7	MAP2K7	O14733
MKK6	MAP2K6	P52564
p38α	MAPK14	Q16539
ERK2	MAPK1	P28482
p90RSK	RPS6KA1	Q15418
CRAF	RAF1	P04049
JNK1	MAPK8	P45983
p110α	PIK3CA	P42336
BRAF	BRAF	P15056
ASK1	MAP3K5	Q99683
p70S6K	RPS6KB1	P23443
MEK1	MAP2K1	Q02750
KRAS	KRAS	P01116
HRAS	HRAS	P01112
SHP2	PTPN11	Q06124
S6	RPS6	P62753

Figure 3

Analysis of dynamics of EGFR signaling. HEK293T cells overexpressing GFP-tagged signaling proteins listed in Table 1 were treated with EGF for 0, 5, 15, 30, and 60 min. Strong abundance-dependent signaling relationships (Supplementary Fig. 10c) are plotted on the signaling map with circle sizes and colors indicating strengths (BP-R2 score) and directionalities (Spearman correlation of bin medians), respectively. The miniaturized network is the same as used in Fig. 2. Overexpression of S6-GFP did not induce any strong signaling relationships (data not shown). For all analyses, data from 3 individual experiment replicates were used.

Figure 4

Analysis of node abundance-dependent EGFR signaling dynamics. (a, b) Schematic plots of amplitude and peak-time analysis. (a) The x-axis (i.e., overexpressed protein as determined by the GFP measurement) was split into 10 bins. (b) Median phosphorylation abundance in each bin was plotted on the y-axis versus time (x-axis) to visualize abundance dependency of signaling dynamics. (c, d) Mass cytometry ion counts (arcsinh transformed, Methods) measured for p-p90RSK (y-axis) as a function of ion counts measured for abundance of CRAF-GFP (x-axis) and EGF stimulation time. The same layouts for (e, f) ERK2-GFP abundance-determined p-p90RSK levels and (g, h) p-ERK1/2 levels are shown. (i) Heat map showing protein abundances with strong influences on signaling amplitudes with color indicating normalized signaling amplitudes. Only overexpressed proteins with an amplitude-ratio higher than 3 fold for more than two of the three replicates were identified as strong influences and are included in the heat map. For (a) to (h), representative examples from the 3 individual experiment replicates are shown. Other replicates are presented in Supplementary File 5. In (i), all replicate data are shown.

In the first analysis, strong and broad signaling responses to overexpression were identified for the upstream kinases PDK1-, GSK3β-, SRC-, and ASK1-GFP without EGF stimulation (Fig. 3. Overall, we identified 59 strong signaling relationships in the unstimulated conditions. Overexpression of many kinases induced strong and positively correlated signaling relationships with their own phosphorylation (Fig. 3, Supplementary File 4). Overexpression of CRAF-, KRAS-, p70S6K-GFP, and others only induced signaling responses upon EGF stimulation (Fig. 3). Notably, under stimulated conditions, KRAS-, CRAF-, and MEK1-GFP levels negatively correlated with phosphorylation levels of downstream kinases p-ERK1/2 and p-p90RSK (Fig. 3). Activating mutations in KRAS and CRAF (Fig. 2), but not protein overexpression alone, may activate oncogenic signaling.

To systematically assess signaling relationships identified by BP-R2, we used the literature-based signaling network, SIGNOR34. For each relationship, we computed the shortest signed directed path length according to the SIGNOR network (Supplementary Table 6). We found that 76% of the strong relationships identified in the unstimulated conditions had paths with a maximum of three steps, highlighting that our approach identifies rather direct signaling relationships. Only 14 abundance-dependent relationships with four or more path steps were identified. Comparison of our strong signaling relationships with literature indicated that many EGF signaling connections that we identified were previously reported. We also propose many relationships that have—to our knowledge—not been previously reported, for example: p90RSK to PDK1 (Ser241), GSK3β to SHP2 (Tyr580), JNK1 to MAPKAPK2 (Thr334), p110α to MKK3 (Ser189), p110α to MKK6 (Ser207), ASK1 to PDK1 (Ser241), ASK1 to GSK3β (Ser9), and ASK1 to AMPKα (Thr172) (Table 2).

Table 2

Relationships with shortest singed directed path length above 3 in the SIGNOR database

Overexpressed POI	Target	Sign	Shortest Signed Directed Path (SIGNOR)	Literature Information
SRC	p-BTK/ITK	1	6	SRC family kinases phosphorylate BTK48
SHP2	p-S6	-1	5	Known regulation49
ASK1	p-PDK1	1	5	Potential novel relationship
SRC	p-PLCγ2	1	5	SRC family kinases activates PLCγ248
ASK1	p-AMPKα	1	4	Potential novel relationship
GSK3β	p-SHP2	1	4	Potential novel relationship
p90RSK	p-PDK1	1	4	Potential novel relationship
JNK1	p-STAT1	1	4	JNK activates STAT150
JNK1	p-MAPKAPK2	1	4	Potential novel relationship
p110α	p-MKK3/6	1	4	Potential novel relationships
HRAS	p-SMAD2/3	1	4	Known crosstalk21
ASK1	p-GSK3β	1	4	Potential novel relationships
PDK1	p-S6	-1	4	Overexpression-induced negative regulation
p70S6K	p-S6	-1	4	Overexpression-induced negative regulation

Phosphorylation levels of many members of the MAPK/ERK pathways showed complex relationships (i.e., measured phosphorylation levels varied over the analyzed POI-GFP range and the relationships did not fit linear, sigmoidal, or quadratic models) with levels of POI-GFPs upon EGF stimulation. These relationships can be explained by abundance-dependent modulation of the signaling dynamics in response to EGF. Thus, in the second set of analyses we examined how signaling dynamics, as quantified by amplitude and peak-time, depended on abundance of an overexpressed protein (Fig. 4). In order to view signaling trajectories as functions of protein abundance, we binned the POI-GFP levels into 10 bins (Fig. 4a, Supplementary File 2). This allowed tracing the signaling trajectories of cells with similar protein overexpression levels (i.e., those in the same bin) over the EGF stimulation time course (Fig. 4b, Supplementary File 5). Strong and robust changes in signaling amplitudes (Fig. 4c-i) and peak-times (Supplementary Fig. 14) were found. Notably, the maximum amplitudes were independent of the overexpression range of a given POI (Supplementary Fig. 15).

We found that high CRAF-GFP and KRAS-GFP abundance strongly reduced signaling amplitudes of p-ERK1/2 and p-p90RSK (Fig. 4c, d, i), whereas high abundance of MEK1-GFP strongly reduced amplitudes and delayed peak-times for p-p90RSK (Fig. 4i, Supplementary Fig. 14). Overexpression of ERK2-GFP led to complex abundance-dependent responses of p-p90RSK and p-ERK1/2 after EGF stimulation (Fig. 4e-h). p-ERK1/2 amplitudes increased and peak-times delayed as a function of ERK2-GFP abundance level (Fig. 4g-i, Supplementary Fig. 14). Intermediate abundance levels of ERK2-GFP also delayed the p-p90RSK peak-times relative to low ERK2-GFP abundance, whereas cells with high ERK2-GFP levels exhibited minimal p-p90RSK signaling dynamics (Fig. 4e, f, i, Supplementary Fig. 14). Overexpression of p90RSK-GFP modulated the signaling amplitude of its potential crosstalk phosphorylation site, Ser241 of PDK1, in an abundance-dependent manner, and increasing expression of p90RSK increased p-PDK1 amplitudes (Fig. 4i). Thus, we observed abundance-dependent signaling dynamics across the range of overexpression levels. Overexpression of upstream signaling proteins (KRAS-, CRAF-, MEK1-, and ERK2-GFP) in the MAPK/ERK pathway led to reduced signaling amplitudes and delayed peak-times of their downstream targets. These observations show that our approach can quantify the role of protein abundance in determining the dynamic signaling response to an extracellular stimulation.

Discussion

Here we present an approach coupling transient overexpression with mass cytometry-based single-cell measurements to characterize signaling network activation states and signaling dynamics over a quasi-continuous, high dynamic range of protein abundance. To highlight the utility of our approach, we present a comprehensive single-cell proteomic analysis of the EGFR network that enabled an analysis of abundance-dependent effects of signaling proteins on state and dynamics of the signaling network. We evaluated the effects of overexpressing 20 EGFR network key nodes with a 60-minute EGF stimulation time course. In each of the 360 conditions, we measured the effect of a POI over a four order of magnitude abundance range on 35 markers by mass cytometry providing a unique and valuable quantitative single-cell resource of abundance dependencies of EGFR signaling.

Previously, the heterogeneity of protein levels after transient transfection was considered problematic. Here, we took advantage of this cell-to-cell variation as it results in a continuous titration of protein abundance over four orders of magnitude. Untransfected cells also provided an internal control for each experiment. We used the multiplexing capabilities of mass cytometry to characterize abundance dependencies of signaling network state and dynamics. Applied to the EGFR signaling network, our approach recapitulated known relationships, suggested previously not described ones, and revealed the intricate modulation of signal amplitudes and peak-times as functions of continuous protein abundance.

Our approach contributes to the understanding of signaling on several levels. First, the approach can be used to study uncharacterized proteins and to suggest additional roles to characterized ones. Second, we were able to directly relate POI abundance with the comprehensive analysis of signaling dynamics in response to stimulation. Such analyses are necessary for understanding of differential signal processing in identical cell types and in disease states characterized by heterogeneity in protein expression such as cancer. Third, the overexpression yields a large dynamic range of signaling activity and can reveal signaling relationships masked by stochastic processes and technical noise under otherwise similar conditions, facilitating the computational analysis of signaling relationships. Fourth, we present a metric termed BP-R2, which allows the quantification of the strengths of arbitrary shaped signaling relationships. BP-R2 was superior to state-of-the-art methods for analysis of our dataset. Fifth, and finally, we were able to infer protein abundance-dependent signaling kinetics from single-cell snapshot data.

Our approach recapitulated known oncogenic signaling behaviors induced by the constitutively active mutants KRASG12V and MEK1DD and identified novel abundance-dependent signaling relationships. For example, p-ERK1/2 was attenuated in cells with highly overexpressed KRASG12V-GFP, potentially due to negative feedback loops or senescence35. Overexpression of the wild-type KRAS-GFP and MEK1-GFP did not induce downstream signaling activation, suggesting that mutations on KRAS or MEK1 are the main drivers of oncogenic signaling. Further, our approach allows study of abundance-dependent signaling dynamics. In the MAPK/ERK pathway, high abundance of upstream signaling mediators KRAS, CRAF, MEK1, or ERK2 reduced amplitudes and delayed peak-times of downstream phosphorylation sites. One possible explanation is that the signal transduction is determined by the competition between active and inactive forms of a signaling protein for substrates. Overexpression increases the total abundance but may reduce the percentage of the active form.

KRAS amplification has been identified in many cancer types. Amplification, however, is not correlated with the phosphorylation of ERK1/236. Rather, KRAS amplification mediates resistance to inhibitors targeting growth pathway related kinases, including EGFR, MET and MEK1/2; KRAS knockdown diminishes the drug resistance37–39. Our results indicate that due to reduced downstream signaling amplitudes in response to EGF stimulation, the dependency of cells on the MAPK/ERK pathway may decrease upon KRAS overexpression, suggesting a mechanism for cancer cell resistance to inhibitors.

Comparing the identified strong signaling relationships with those in the SIGNOR database, we propose previously not described signaling relationships, e.g.: 1) Our data suggest that p90RSK potentially forms a positive feedback loop and activates the upstream signaling protein PDK1. 2) GSK3β has been identified as a central signaling controller and has multiple substrates40; our results suggest that SHP2 is a potential direct or indirect target of GSK3β. 3) We also propose that JNK1 is a MAPKAPK2 activator. 4) PI3K and MKK3/6 are known to be regulated by RAC141; our results suggest PI3K activates MKK3/6 independently. 5) Recent studies indicate that ASK1 contributes in negative regulation of PDK1 through phosphorylation on Thr254 of PDK142; We observe ASK1 overexpression-induced PDK1 phosphorylation on Ser241, inducing PDK1 activity and downstream GSK3β phosphorylation on Ser9. 6) In addition to the known AMPK-mediated ASK1 activation43, our data indicates ASK1 activation of AMPKα via phosphorylation on Thr172. 7) We have also observed negative correlations between the abundance of p70S6K or PDK1 to the phosphorylation level of S6 (Ser235/Ser236), indicating overexpression-induced-negative feedback regulations.

Our method has several limitations. First, we do not measure the endogenous expression level of the POI. However, exogenous expression is linearly correlated with the total protein level (Supplementary Fig. 4a), validating GFP as readout of the total POI. Second, all results in mass cytometry rely on antibodies; for this work, all antibodies were thoroughly validated (Supplementary Table 3). Third, we do not measure the abundance range of the studied proteins in cancer cells, however, proteome studies of cancer cells and databases such as PaxDb11 indicate a range similar to those studied here. Fourth, high expression levels of a protein kinase may induce non-specific phosphorylation; however, our data allows choosing the analyzed expression range in silico, thus such effects can be excluded.

The approach described here provides a method to study how the abundance variance of signaling proteins in different tissues and cell lines results in distinct signaling behaviors. The application of our approach to synthetic biology, stem cell biology, developmental biology, and cancer-related processes, such as the epithelial-mesenchymal transition, will enable quantitative identification of key proteins and signaling determinants in cell differentiation at phenotypical switching points. We envision that determining which signaling relationships and thresholds enable diseased cells to overcome drug treatment will be a highly relevant application.