Receptor-based mechanism of relative sensing and cell memory in mammalian signaling networks

Detecting relative rather than absolute changes in extracellular signals enables cells to make decisions in constantly fluctuating environments. It is currently not well understood how mammalian signaling networks store the memories of past stimuli and subsequently use them to compute relative signals, that is perform fold change detection. Using the growth factor-activated PI3K-Akt signaling pathway, we develop here computational and analytical models, and experimentally validate a novel non-transcriptional mechanism of relative sensing in mammalian cells. This mechanism relies on a new form of cellular memory, where cells effectively encode past stimulation levels in the abundance of cognate receptors on the cell surface. The surface receptor abundance is regulated by background signal-dependent receptor endocytosis and down-regulation. We show the robustness and specificity of relative sensing for two physiologically important ligands, epidermal growth factor (EGF) and hepatocyte growth factor (HGF), and across wide ranges of background stimuli. Our results suggest that similar mechanisms of cell memory and fold change detection may be important in diverse signaling cascades and multiple biological contexts.


Introduction
In biological systems, concentrations of extracellular signaling molecules, such as hormones and growth factors, often vary by orders of magnitude. Therefore, the ability to sense relative rather than absolute signals, that is detect fold changes in extracellular cues, is critical for making accurate decisions in different biological contexts (Alon, 2019). Relative sensing requires both the ability to store memories of past environmental stimuli and the capacity to quickly and efficiently compute relative signals (Adler and Alon, 2018).
Relative sensing of environmental inputs has been previosuly investigated in bacteria, with the E. coli chemotaxis being a classic example (Mesibov et al., 1973;Barkai and Leibler, 1997;Alon et al., 1999;Shoval et al., 2010). Studies have also explored relative sensing in a variety of eukaryotic systems. When responding to constant stimuli, experiments with the signaling proteins ERK (Cohen-Saidon et al., 2009) and b-catenin (Goentoro and Kirschner, 2009) showed that fold changes in their nuclear activity were robust to cell-to-cell variability (Cohen-Saidon et al., 2009) and variability in signaling network parameters (Goentoro and Kirschner, 2009). These observations suggested that gene expression of target genes may respond, at the single cell level, to fold changes rather than absolute activities of these proteins. Later studies of the NF-kB (Lee et al., 2014) and TGF-b/SMAD pathways (Frick et al., 2017) also showed that genes directly controlled by these proteins often respond to their fold changes at the single cell level. Recent work has explored relative sensing at the organism level in plants, where the chlorophyll activity was found to be proportional to the fold change in external light intensity (Tendler et al., 2018).
Despite the insights gained in the aforementioned studies, the molecular mechanisms allowing cells to detect fold changes in extracellular stimuli are not well understood. The key unresolved questions are: (1) where and how the memories of background extracellular stimuli are stored within the cell, (2) what makes these memories specific to particular stimuli, and (3) how the cells subsequently use the stored memories to compute fold changes.
In this work, using the growth factor-activated PI3K/Akt signaling pathway, we describe a novel non-transcriptional mechanism of relative sensing in mammalian cells. The mechanism operates on fast timescales of dozens minutes to hours, and across more than an order of magnitude of extracellular background stimuli. We derive key aggregate parameters of the signaling cascade that determine the accuracy and the background range of relative sensing. We also experimentally validate the accuracy of relative sensing by stimulating cells with multiple fold changes of two physiologically important ligands, EGF and HGF. Furthermore, we demonstrate that ligand relative sensing is reliably propagated to an important downstream target of the PI3K/Akt pathway.
To understand how the immediate-early dynamics of the Akt pathway depend on the background level of growth factors, we used immunofluorescence to quantify the levels of pAkt in epidermal growth factor (EGF)-stimulated human non-transformed mammary epithelial MCF10A cells (Materials and methods, Figure 1-figure supplement 1). Within minutes of continuous stimulation with EGF pAkt reached maximum response, and then decayed to low steady state levels within hours ( Figure 1a). The resulting steady state pAkt levels were approximately independent of the EGF stimulus, indicating an approximately adaptive response (Friedlander and Brenner, 2009;Shoval et al., 2010; Figure 1-figure supplement 2). In the sensitive range of EGF concentrations, maximal pAkt response was approximately proportional to the logarithm of the EGF stimulus ( Figure 1b). Quantitative western blot experiments demonstrated that in this logarithmic regime, pAkt levels were approximately linearly proportional to the phosphorylation level of EGF receptors (EGFRs) (Figure 1-figure supplement 3). The logarithmic dependence of EGFR phosphorylation levels on EGF stimulation has been previously attributed to a mixture of receptor species with varying affinities to the ligand, negative cooperativity of ligand binding to receptor dimers, and oligomeric aggregation of receptors (Kawamoto et al., 1983;Chatelier et al., 1986;Wofsy et al., 1992;Macdonald and Pike, 2008;Huang et al., 2016).
Continuous stimulation with EGF resulted in the abundance of cell-surface EGF receptors (sEGFR) also decreasing proportionally to the logarithm of the background EGF level, and reaching a new steady state within hours ( Figure 1c). Notably, prior exposure with EGF desensitized cells to subsequent EGF stimulations in a quantitative manner. When we first pre-exposed cells to different levels of EGF for 3 hr and then stimulated them with the same final EGF concentration (2 ng/ml)), the maximal pAkt response decreased monotonically with increasing pre-exposure EGF levels ( Figure 1d). These experiments demonstrate that the pAkt response to an abrupt EGF stimulation is strongly affected by background EGF levels, and that this effect is likely mediated by the endocytosisbased removal of activated EGFRs from the cell surface (Wiley et al., 1991).
Using pharmacological perturbations of the EGFR/Akt pathway, we confirmed that the desensitization of the phosphorylation response ( Figure 1d) was likely due to receptor-based mechanisms upstream of Akt activation, and did not depend on its downregulation, for example, through phosphorylation-dependent Akt degradation (Wu et al., 2011). Specifically, we used SC79, a small molecule which promotes Akt phosphorylation even in the absence of extracellular ligands (Jo et al., 2012). Unlike the desensitization observed in the growth factor-induced pAkt response (Figure 1d), the pAkt response following stimulation with SC79 did not depend on the background EGF preexposure ( Figure 1-figure supplement 4). This result supports the conclusion that the EGF desensitization mechanism was upstream of Akt.
To understand how background EGF levels affect the pAkt response to subsequent EGF stimulation we next constructed an ordinary differential equation (ODE) model of EGF-dependent Akt phosphorylation. The model included several well-established features of the EGFR signaling cascade (Chen et al., 2009), such as endocytosis and degradation of activated receptors (Materials and methods) (Figure 2a). We constrained the ranges of model parameters based on literaturederived estimates (Supplementary file 1a), and fitted the model using experimental data on pAkt time courses (Figure 1a) and steady state sEGFR levels ( Figure 1c) at different doses of EGF (d) Desensitization of the maximal pAkt response to an abrupt EGF stimulation. MCF10A cells were pre-treated with increasing background doses of EGF (x-axis) for three hours, followed by a second abrupt stimulation with the same concentration of EGF (2 ng/ml); the inset shows a schematic illustration of the experimental protocol. In all subpanels, error bars represent the standard deviation of n = 3 technical replicates. Source data: pakt_timecourses_first_step.mat and segfr_150_180mins.doseresponse. mat (available in Source code 1). The online version of this article includes the following figure supplement(s) for figure 1:    stimulations. We then used simulated annealing to optimize model parameters (Materials and methods, Figure 2-figure supplement 1), and considered multiple distinct parameter sets from the optimization runs for further computational analysis.
Using the fitted dynamical model (Figure 2-figure supplement 1), we explored the ability of the Akt pathway to respond to relative, rather than absolute, changes in EGF levels. To that end, we simulated the pAkt response by exposing the model in silico to a range of background EGF levels followed by different abrupt fold change increases in EGF concentration (Figure 2b). The model predicted that the maximal pAkt response indeed depends primarily on the EGF fold change relative to Schematic of the computational model of the EGFR signaling cascade leading to phosphorylation of Akt. Rate constants marked with asterisks correspond to reactions associated with activated (phosphorylated) receptors. Only a subset of reactions in the network are shown for brevity. (b) In silico protocol used to explore relative sensing, showing the temporal profiles of EGF stimulation (top) and the corresponding profiles pAkt response (bottom). Cells were first exposed to various background EGF stimulations (blue and red) and were next subjected to the same abrupt fold change in EGF at time t 0 . The resulting maximal pAkt responses were similar for the same EGF fold change independent of background EGF stimulation, indicating relative sensing. (c) The maximal pAkt response observed after exposing the ODE model in silico to different background EGF levels (x axis), followed by a 2-, 3-, 4-, or 6-fold increase (different colors) in EGF; inset shows pAkt response over a wider range of background EGF levels. (d) Maximal pAkt responses (y axis) induced by stimulation with different EGF background levels (data points with the same shape and color) were combined and plotted as a function of the EGF fold change (x axis). Dashed line represents log-linear fit to data (Pearson's r 2 = 0.96, regression p value < 10 À15 ). In all subpanels, error bars represent the standard deviation of n = 10 model fits. Source code: https://github.com/ dixitpd/FoldChange/. The online version of this article includes the following figure supplement(s) for figure 2: the background stimulation levels ( Figure 2c). This relative sensing of EGF stimuli occurred over an order of magnitude of background EGF concentrations, and the resulting pAkt response was approximately proportional to the logarithm of the EGF fold change (Figure 2d). Notably, the model predicted relative sensing exactly in the range of EGF background concentrations where sEGFR endocytosis was sensitive to background ligand stimulation. At low EGF background concentrations (<0.01 ng/ml), no substantial sEGFR removal was predicted at the steady state (Figure 2figure supplement 1), and consequently there was no significant desensitization of the pAkt response. In that regime, the pAkt response to an abrupt fold change depended primarily on the absolute EGF level. In contrast, at high background EGF concentrations (>1 ng/ml), a large fraction of sEGFR was already removed from cell surface and consequently the network responded only weakly to further EGF stimulation.
Next, we experimentally tested the model-predicted relative sensing in MCF10A cells. Cells were first treated with various background EGF concentrations for three hours to ensure that sEGFR reached steady state levels (Figure 1c), and that pAkt had decayed after a transient increase ( Figure 1a). As in the computational analysis (Figure 2b), cells were then exposed to different fold changes in EGF levels; pAkt levels were measured at 2.5, 5, 10, 15, 30 and 45 min after the step increase in EGF stimulation (  The maximal pAkt responses after exposing cells to different background EGF levels (x axis) for 3 hr, followed by 2-, 3-, 4-, and 6-fold increases (different colors) in EGF. Inset shows experimental pAkt response over a wider range of background EGF levels. (b) Maximal pAkt responses (y axis) to fold changes in EGF depended approximately logarithmically on the fold change. Maximal pAkt responses induced by stimulation with various EGF background levels (data points with the same shape and color) were combined and plotted as a function of the EGF fold change (x axis). Dashed line represents log-linear fit to the data (Pearson's r 2 = 0.93, regression p value < 10 À11 ). In all subpanels, error bars represent the standard deviation of n = 3 technical replicates. Source data: expt_data.mat (available in Source code 1). The online version of this article includes the following figure supplement(s) for figure 3:       pAkt responses. The concentration range in which we obsered relative sensing was consistent with recent estimations of in vivo EGF levels (Pinilla-Macua et al., 2017). In close agreement with the computational model predictions, the maximal pAkt response was approximately proportional to the logarithm of EGF fold change ( Figure 3b). Interestingly, in addition to the maximal pAkt response, approximate relative sensing was also observed for the time integral of pAkt levels ( To better understand the mechanism responsible for the observed relative sensing of extracellular EGF concentration, we next constructed a simplified analytical model of the signaling network (see Appendix). This model revealed that, across a broad range of background concentrations, the steady-state abundance of cell surface receptors R ½ T decreases approximately log-linearly as a function of the background ligand (EGF) concentration L ½ 0 (Equation 1 and Figure 4a): and that the maximal receptor phosphorylation response LR Ã 2 Â Ã max depends approximately log-linearly on the level of the subsequent stimulation L ½ 1 and linearly on the steady-state receptor abundance R ½ T (Equation 2, Figure 4b): where a and b are numerical constants (Appendix). As a result of these relationships, the phosphorylation response LR Ã 2 Â Ã max after an increase of ligand concentration from L ½ 0 to L ½ 1 depends, in agreement with computational and experimental analyses, approximately on the logarithm of the stimulation fold change The analytical model (Appendix) also revealed that the range of the background ligand concentrations where relative sensing is observed is primarily determined by two aggregate systems parameters, which we denote a and b (Equations 4 and 5). The parameter a quantifies the ability of the signaling network to capture the input signal (EGF) and elicit a downstream phosphorylation response. The parameter b quantifies the ability of the network to preferentially internalize and degrade active (phosphorylated) receptors relative to inactive (non-phosphorylated) receptors. The two aggregate parameters are expressed as follows: where k p is the rate of receptor phosphorylation and k dp is the rate of receptor de-phosphorylation, k 2 is the rate of receptor dimerization, k À2 is the dissociation rate of receptor dimers, and R 0 is the total number of cell-surface receptors at the steady state in the absence of extracellular stimuli and and k i ; k rec ; k deg are correspondingly the rates of internalization, recycling, and degradation of the active (phosphorylated) and non-active receptors. Notably, an increase in the value of a increases signal sensitivity and receptor dimerization and phosphorylation. This shifts the relative sensing range to lower ligand concentrations ( Figure 5a). In turn, an increase in the value of b increases the fraction of active receptors being internalized and degraded. This increases the range of background ligand concentrations where the relative sensing is observed ( Figure 5b). Based on the best-fit ODE model parameter sets, we estimate a~15 and b~40 (Appendix). As an example, in Figure 5 we show the scaled phosphorylation response to a six-fold change in EGF concentration as a function of the scaled background ligand concentration u 0 for different values of a ( Figure 5a) and b ( Figure 5b); the green arrows in the figure represent the predicted range of fold change detection. The model analysis showed that the relative sensing occurs across over an order of magnitude of background ligand concentrations (Appendix). Furthermore, the analytical model revealed that relative sensing does not require receptor dimerization, and similar sensing mechanisms can operate in pathways where signaling is initiated by monomeric receptors (Appendix).
In addition to EGF, Akt phosphorylation can be induced by multiple other ligands, including hepatocyte growth factor (HGF) (Stuart et al., 2000) which binds to its cognate receptor cMet Figure 5. Analytical model predicts the range of approximate EGF relative sensing. Scaled phosphorylation response to a six-fold change in extracellular EGF concentration as a function of the scaled background ligand concentration (u 0 ¼ ½L 0 =K d1 ). (a) The phosphorylation response as a function of background ligand concentration (x-axis) is shown for different values of the parameter a (y-axis), when the parameter b is fixed at b=40. (b) The phosphorylation response as a function of background ligand concentration (x-axis) is shown for different values of b (y-axis), when a is fixed at a=15. The colors represent the scaled phosphorylation response. The green dashed lines delineate the parameter ranges where the fold change detection is >90% accurate. The horizontal black lines correspond to the parameter values a=15 and b=40, which were estimated from experimental data fits; the horizontal green double arrows represent the predicted range of relative sensing for the investigated PI3K-Akt cascade. (Viticchiè and Muller, 2015). Similar to EGFRs, upon ligand binding, cMet receptors dimerize (Kong-Beltran et al., 2004) and cross-phosphorylate each other; this leads to phosphorylation of multiple downstream targets, including Akt. To investigate the specificity of the receptor-based cell memory to past ligand exposures, we used the two ligands, EGF and HGF, which share many signaling components downstream of their cognate receptors (Xu and Huang, 2010). We exposed cells to background doses of either HGF or EGF for three hours, and then stimulated cells using either the same or the other growth factor to elicit pAkt response (Figure 6a,b). Pre-exposure with HGF did not substantially downregulate EGF-induced pAkt responses, but substantially Figure 6. Desensitization and ligand-specific cell memory for EGF-and HGF-induced pAkt responses. MCF10A cells were first exposed to various background concentrations of either HGF or EGF for three hours, and then abruptly stimulated using either the same or the other growth factor. pAkt levels were then measured 10 min after the addition of the second stimulus. (a) EGF-(blue, 2.5 ng/ml) or HGF-(red, 4 ng/ml) induced pAkt response in cells pre-exposed with various background concentrations of HGF (x axis). (b) EGF-(blue, 2.5 ng/ml) or HGF-(red, 4 ng/ml) induced pAkt response in cells pre-exposed with various background concentrations of EGF (x axis). (c) The maximal pAkt response in MCF10A cells exposed to different background doses of HGF (x axis) for 3 hr, followed by 2-, 4-, and 8-fold increase (different colors) of HGF. Inset shows experimental pAkt response over a wider range of background HGF levels. (d) The maximal pAkt responses (y axis) to HGF fold changes depended approximately logarithmically on the fold change (x axis). Maximal pAkt responses induced by stimulation with various HGF background levels (data points with the same shape and color) were combined and plotted as a function of the HGF fold change (x axis). Dashed line represents log-linear fit to data (Pearson's r 2 = 0.88, regression p value < 10 À6 ). In all subpanels, error bars represent the standard deviation of n = 3 technical replicates. Source data: expt_data.mat (available in Source code 1). The online version of this article includes the following figure supplement(s) for figure 6:  decreased HGF-induced responses (Figure 6a). Similarly, we observed a relatively small desensitization of HGF-induced responses due to pre-exposure with EGF, while there was a significant desensitization of EGF-induced pAkt responses (Figure 6b). We further confirmed that exposure of MCF10A cells to various concentrations of HGF leads to pronounced HGF-dependent removal of cMet from the cell surface, without significant removal of sEGFR ( Figure 6-figure supplement 1). Similarly, the pre-exposure of cells to EGF leads to EGF-dependent removal of sEGFR without a significant change in surface cMet abundance ( Figure 6-figure supplement 1). These observations support the mechanism in which the relative sensing of extracellular ligands relies on the memory of their past exposures effectively encoded in the abundances of their cognate cell-surface receptors.
Given the observed HGF-dependent removal of cell surface cMet receptors and the resulting pAkt desensitization, we investigated next whether the maximal pAkt response depends, similarly to EGF, on the relative fold changes in the level of extracellular HGF. To that end, we exposed cells to a range of different background levels of HGF, and then stimulated cells with different fold changes in HGF concentrations (Figure 6c,d and Figure 6-figure supplement 2). These experiments demonstrated that HGF-induced phosphorylation of Akt also depends primarily on the fold change in extracellular HGF concentration across almost an order of magnitude of background HGF exposures (between 0.1 and 1 ng/ml HGF) ( Figure 6c). Moreover, like EGF, the maximum pAkt levels depended approximately log-linearly on the HGF fold change (Figure 6d).
Relative sensing of extracellular ligands should affect important downstream biological targets of the PI3K-Akt pathway. The FoxO3 transcription factor is a key effector of the pathway, and it is involved in diverse cellular processes including apoptosis, proliferation, and metabolism (Webb and Brunet, 2014). Akt phosphorylation of FoxO3 leads to its translocation from the nucleus to cytoplasm and subsequent transcriptional deactivation (Webb and Brunet, 2014). Notably, following Akt activation, the typical nuclear translocation timescale for FoxO family proteins is short (less than 5 min) (Gross and Rotwein, 2017). To investigate FoxO3 activation induced by EGF stimulation, we used quantitative immunofluorescence to measure its nuclear-to-cytoplasm ratio (Worster et al., 2012). We exposed cells to two different background EGF levels for three hours, and then treated them with two different abrupt fold changes in EGF concentrations. Consistent with relative sensing by pAkt, the nuclear-to-cytoplasmic ratio of FoxO3 also reflected the relative, rather than the absolute changes in EGF stimulation (Figure 7 and Figure 7-figure supplement 1). Thus, relative sensing of the growth factor signal is faithfully transmitted in MCF10A cells to at least some of the physiologically important effectors of the PI3K-Akt pathway.

Discussion
Receptor endocytosis and down-regulation, following ligand stimulation, has been canonically associated with signal and circuit desensitization (Friedlander and Brenner, 2009;Sorkin and von Zastrow, 2009;Ferrell, 2016). Our study suggests an additional and more quantitative role of receptors endocytosis in mammalian cells. Specifically, receptor endocytosis may allow cells to continuously monitor signals in their environment (Becker et al., 2010;Brennan et al., 2012;Mitchell et al., 2015) by dynamically adjusting the number of ligand-cognate receptors on the cell surface. Our analysis also demonstrates that the memory of past stimuli, effectively encoded in the number of surface receptors, may be signal-specific, at least for some ligands, due to the selective removal of ligand-cognate receptors. The combination of logarithmic pAkt response, and the logarithmic dependence of the ligand-specific memory on the background signal, allows cells to respond to relative changes in environmental stimuli. We note that the described relative sensing mechanism is not a direct consequence of either simple adaptation to various levels of background signals or logarithmic activation response (Shoval et al., 2010;Adler et al., 2017;Adler and Alon, 2018).
Previous studies (Cohen-Saidon et al., 2009;Goentoro and Kirschner, 2009;Lee et al., 2014) have demonstrated that transcriptional motifs may efficiently buffer cell-to-cell variability in signaling components when responding to a constant extracellular stimulation. In contrast, our study describes a non-transcriptional mechanism of sensing extracellular signal changes relative to past extracellular stimulation. Although the pAkt response to an abrupt stimulation is relatively fast (~5-15 min, Figure 1a), and therefore non-transcriptional in nature, the sustained production and delivery of cell surface receptors is essential to establishing the signal-dependent and receptor-mediated memory. Therefore, sustained transcription and translation of network comonents are necessary for proper functioning of the described sensing mechanism.
Although there are usually~10 5 -10 6 EGFR receptors on mammalian cell surface , the downstream network response, for example Akt phosphorylation, often saturates when only a relatively small fraction (5-10%) of the receptors are bound to their cognate ligands (Chen et al., 2009;Shi et al., 2016). Our study suggests that one potential advantage of such a system architecture is that, beyond simple signal activation, it may endow cells with a large dynamic range of receptor abundances to memorize stimulation levels of multiple extracellular ligands (Hart et al., 2013;Nandagopal et al., 2018). Notably, signal-mediated removal has been reported for many other receptors and signaling systems, such as the G protein-coupled receptors (GPCRs) (Ferguson, 2001), involved in various sensory systems, and AMPA-type glutamate receptors (Guskjolen, 2016), implicated in synaptic plasticity. Therefore, similar relative sensing mechanisms may be important in multiple other receptor-based signaling cascades and across different biological contexts.

Measurement of EGF signaling responses
MCF10A cells were obtained from the ATCC and grown according to ATCC recommendations. Cell identity was confirmed by short tandem repeat (STR) profiling at the Dana-Farber Cancer Institute and cells were tested with the MycoAlert PLUS mycoplasma detection kit (Lonza) and found to be free of Mycoplasma prior to analyses. For experiments, 96 well plates (Thermo Scientific) were coated with type I collagen from rat tail (Sigma-Aldrich) by incubating plates with 65 ml of 4 mg/ml collagen I solution in PBS for 2 hr at room temperature, washed twice with PBS using a EL406 Figure 7. Relative sensing of EGF concentrations by pAkt is propagated to FoxO3. MCF10A cells were first exposed to two background concentrations of EGF for three hours, and then were stimulated with 3-and 6fold increase in EGF concentrations. The ratio of nuclear-to-cytoplasmic FoxO3 levels (y-axis) was measured using quantitative immunofluorescence (Materials and methods) after 15 min of the EGF fold changes. Statistical significance was calculated using the Wilcoxon rank sum test (n = 5); * corresponds to p<0.01, and n. s. corresponds to p>0.1. Error bars represent the standard deviation of n = 5 technical replicates. Source data: expt_data.mat (available in Source code 1). The online version of this article includes the following figure supplement(s) for figure 7: Microplate Washer Dispenser (BioTek), and then sterilized under UV light for 20 min prior to use. Cells were harvested during logarithmic growth and plated into collagen-coated 96 well plates using a EL406 Microplate Washer Dispenser. Cells were grown in 200 ml of complete medium for 24 hr, serum starved twice in starvation media (DMEM/F12 supplemented with 1% penicillin-streptomycin and 0.1% bovine serum albumin), incubated in 200 ml of starvation media for 19 hr, washed twice more, and incubated in 200 ml of starvation media for another hour. This time point constituted t = 0 for all experiments. Treatment solutions were created by manual pipetting or by dispensing the appropriate amounts of epidermal growth factor (EGF, Peprotech), hepatocyte growth factor (HGF, Peprotech), or SC-79 (Sigma) into starvation media using a D300 Digital Dispenser (Hewlett-Packard). At t = 0 cells were stimulated with 100 ml of treatment solution and then incubated for the indicated times. For experiments requiring a second stimulus, cells were treated with an additional 100 ml of treatment solution at 3 hr and incubated for the indicated times. For fixation, 100 ml of supernatant were removed from the wells, replaced by 100 ml of 12% formaldehyde solution (Sigma) in phosphate buffered saline (PBS), and incubated for 30 min at room temperature.
All subsequent washes and treatments were performed with the EL406 Microplate Washer Dispenser. Cells were washed twice in PBS and permeabilized with 0.3% Triton X-100 (Sigma-Aldrich) in PBS for 30 min at room temperature (this step was omitted for measuring the surface expression of cMET and EGFR), washed once again in PBS, and blocked in 40 ml of Odyssey blocking buffer (LI-COR Biotechnology) for 60 min at room temperature. Cells were incubated with 30 ml of anti-phospho-Akt (Cell Signaling Technologies, 4060, 1:400), FoxO3 (Cell Signaling Technologies, 2497, 1:200), anti-Met (R and D Systems, AF276, 1:150), or anti-EGFR (Thermo Fisher Scientific, MA5-13319, 1:100) over night at 4˚C. Cells were washed once in PBS and three times in PBS with 0.1% Tween 20 (Sigma-Aldrich; PBS-T for 5 min each and incubated with 30 ml of a 1:1000 dilution of secondary antibodies conjugated with Alexa Fluor 647 in Odyssey blocking buffer for 60 min at room temperature. Cells were washed two times in PBS-T, once with PBS, and stained for 30 min at room temperature with whole cell stain green (Thermo Fisher Scientific) and Hoechst (Thermo Fisher Scientific). Cells were washed three times in PBS, covered in 200 ml of PBS, and sealed for microscopy. Cells were imaged using an Operetta microscope (Perkin Elmer).

Image Processing
Images were analyzed using the Columbus image data storage and analysis system (Perkin Elmer) to quantify single cell fluorescence measurements from each imaged well. The reported intensity values were obtained by first subtracting the background fluorescence of the well and subsequently the levels of pAKT at no stimulation at the same time. From each well we thus obtained a distribution of single cell measurements of a given target (pAkt, FoxO3, scMET or sEGFR). In each distribution we discarded the top and bottom 5% of points to remove outliers due to imaging and detection errors. The nuclear FoxO3 to cytoplasmic FoxO3 compartmentalization ratio was determined by the mean intensity in each area after image segmentation based on Hoechst and whole cell stain green at the single cell level. After that, we calculated the average of the resulting single cell distributions. For each condition, we performed multiple technical repeats (multiple wells), and as a final result reported the average of the corresponding single cell distribution averages and the associated standard deviations.

Computational methods
In this section, we describe in detail (1) the model of the EGF/EGFR/Akt signaling pathway, (2) model assumptions, (3) model parameters and their bounds, (4) various relevant biological constraints that were imposed while fitting the model to the data, (5) the error function that was minimized in the parameter search, (6) the numerical procedure used to minimize the error function between model predictions and experimental data, and (7) in silico predictions.

General structure of the computational model
The dynamic ODE model describing the EGF/EGFR signaling cascade leading to Akt phosphorylation (Supplementary file 1a, b, and Equations A2-A20) was based on the previous work by Chen et al. (2009). We retained the components of the model relevant to EGF-dependent phosphorylation of EGFR and the subsequent cascade responsible for Akt phosphorylation. The resulting model consisted of 20 chemical species (see Supplementary file 1a) and was described by 24 parameters (20 reaction rate constants and four total species concentrations, Supplementary files 1-3, and Equations A2-A20).
The model included processes across three cellular compartments: cell surface (plasma membrane), cytoplasm, and endosomes. The model included interactions of the ligand with the receptors (ligand-binding and unbinding to receptor monomers and dimers) and subsequent receptor dimerization and undimerization. The model also included internalization of phosphorylated and unphosphorylated receptors, their recycling and degradation, phosphorylation and dephosphorylation by phosphatases.

Main assumptions of the model
In agreement with available literature (Wiley and Cunningham, 1982;Herbst et al., 1994), we assumed that the rates of internalization, recycling, and degradation are different for inactive (unphosphorylated) and active (phosphorylated) receptors (Supplementary file 1a). We assumed that EGFR phosphatases in MCF10A cells are present at exceedingly high concentrations (Kleiman et al., 2011), and therefore we implemented the corresponding reaction of dephosphorylation of phosphorylated EGFRs (pEGFRs) as a first order reaction. We further assumed that activated receptors on plasma membrane and in endosomes are dephosphorylated by the phosphatase with the same rate (Kleiman et al., 2011).
We implemented PIP2 phosphorylation by pEGFR on the plasma membrane as a simplified effective first order process. Following the receptor-driven phosphorylation of PIP2 we retained the canonical signaling cascade of the PI3K/Akt activation ( Figure 2a in the main text). We also implemented a first order reaction for action of the phosphatase on pAkt.
We assumed that cells are at steady state in terms of the abundances of ligand-free cell surface and endosomal receptors prior to ligand exposure. Specifically, prior to ligand exposure, the number of ligand-free EGFR monomers on cell surface and in endosomes, were derived based on the steady state condition of the corresponding equations.
In agreement with the literature (Haugh and Meyer, 2002;Park et al., 2003), we assumed that Akt can be phosphorylated only by cell-surface pEGFR, and not by endosomal pEGFR. Finally, we assumed that over the course of simulation extracellular ligand concentration remained constant, unless a step increase in EGF was applied. Here we refer to the background ligand stimulation as stimulation applied at time t = 0 to the cells that were previously not exposed to the ligand.

Model parameters
The model parameters consisted of 4 total species abundances (PIP2, Akt, PDK1, and EGF receptors) and 20 rate constants. We collected multiple values of these parameters from literature (Supplementary file 1a). In our search for optimal rate parameters fitted to data, we allowed rate parameters to vary within half an order of magnitude from the lowest and the highest literature derived estimate (Supplementary file 1a). For parameters, for which experimental estimates were not available, we allowed up to four orders of magnitude in variation.
In addition, we allowed one and a half orders of variation in first order rate of EGF unbinding from receptors and the rate of EGFR phosphorylation in order to account for spatial organization of the receptors on the cell surface (Mayawala et al., 2006). We fixed the rate of pEGFR phosphatase according to the measurement of this constant in MCF10A cells (Kleiman et al., 2011). In accordance with literature parameter estimates (Supplementary file 1a), we constrained the rate of ligand unbinding, receptor undimerization, receptor phosphorylation, and receptor dephosphorylation to be at least 10 times faster than receptor internalization (Wiley et al., 1991;Herbst et al., 1994;Chen et al., 2009;Kleiman et al., 2011).

Additional constraints
In addition to the constraints imposed on network parameters directly through the experimentally measured data at EGF stimulations, we also required several additional constraints to better capture biology of EGFR signaling based on known literature. These constraints were either added as 'hard' constraints: parameter sets that did not agree with hard constraints were rejected, or as 'soft' constraints: parameter sets that did not agree with soft constraints were penalized using additional terms in the error function.

Hard constraints
We constrained the total number of EGF receptors prior to EGF exposure to be between 10 5 -10 6 per cell, and surface EGFR to be within 10 5 -10 6 receptors per cell, in agreement with EGFR abundances reported for MCF10A cell lines (Niepel et al., 2013). In the model, the number of cell surface receptors was not a free parameter, but was calculated based on the steady state condition of the differential equations that describe the system.

Soft constraint
We also implemented a 'soft' constraint that ensured realistic levels of phosphorylated Akt molecules. We required that at least 10% of total Akt gets phosphorylated at EGF doses close to Akt saturation, that is in our case 3.16 ng/ml EGF (Chen et al., 2009).

Error function
The error function quantified the disagreement between model predictions and data and the soft constraints. At a given point in parameter space, we solved the system of ODEs describing EGFinduced Akt phosphorylation (see below) using the MATLAB and obtained model solutions {S i } for all experimentally measured conditions. Importantly, while the model predicts protein concentrations in units of number of molecules per cell, our experiments measured protein concentration up to a scaling factor. Therefore, we rescaled the model prediction using maximum likelihood linearregression estimate (MLE) for both pAkt and sEGFR data between the model and the data respectively. Specifically, separately for pAkt and sEGFR, we fitted a linear model between the predictions from the ordinary differential equation model and the corresponding experimental measurements across multiple EGF doses and time points. We rescaled the predictions based on the slope and the intercept of the linear model fit. The scaled predictions were used in the evaluation of the error.
The error function comprised of two different contributions. The first term was defined as the sum of the squared differences between the model predictions {S i } at parameter value and the corresponding experimental data D, taking into account corresponding experimental errors s s (Equation S1a). Next, we imposed the soft constraints described above as squared error terms (Equation S1b, and Supplementary file 1b for species abbreviation). The total error function was the sum of these two contributions (Equation S1c).
The standard deviation 0.004 in Equation S1b was chosen to ensure that the maximum pAkt levels were guaranteed to be above 10% of total Akt levels. Lower values lead to a very high rejection rate in the simulated annealing procedure and higher values were likely to return parameter points that did not satisfy the constraint that maximum pAkt levels were at least 10% of total Akt levels. The active endocytosis and degradation of cell surface receptors in our system occurred mostly between EGF doses of 0 ng/ml and 3.16 ng/ml. Accordingly, we fit the model using experimental data collected in the same range of EGF stimulations.
In Equation 1 a, index k runs through all n experimentally measured data points (5 doses x 7 time points = 35 total points) and sEGFR measurements (6 doses x 2 time points = 12 total points).
Overall, the error function had a total of 50 terms (35 pAkt measurements, 12 sEGFR measurements, and one soft constraint). We minimized this error by searching through the parameter space using simulated annealing (SA) described in the next section.

Simulated Annealing optimization
Given that the mechanistic ODE models constrained by experimental measurements of several dynamical quantities are usually underdetermined (Chen et al., 2009), we used simulated annealing (SA) (Kirkpatrick et al., 1983), to numerically search the model's parameter space.
The overall error (Equation S1c) was minimized with SA in order to determine parameter sets that are most consistent with the experimental measurements. Following standard SA optimization scheme, we ran a random walk in the model's parameter space. At each point in the parameter space we accept or reject a next proposed parameter set according to the Metropolis criterion and a likelihood that is the negative exponential of the error function in Equations S1a. Following a conventional SA protocol, we used an additional parameter, temperature, which allowed steps with relatively large change in the likelihood score to explore large parameter space. The temperature was decreased gradually to find a local minimum of the likelihood function.
To find multiple parameter sets that fit the experimental data, we ran 100 independent SA chains with randomly selected starting points spread out across allowed parameter ranges. Each chain was started at high temperature and was cooled down in 12 stages to the lowest temperature (using the sequence of temperatures: 400, 200, 100, 50, 20, 10, 5, 2, 1, 0.5, 0.25, 0.1). At each temperature, 1500 steps in parameter space were performed. In each step, on an average four randomly chosen parameters (out of the 24) were changed in order to speed up the search in the parameter space.

Predictions from SA
For individual chains, the parameter set with the lowest error was recorded. The averages of parameter values from the 10 best-fit chains are shown in Supplementary file 1a. We used these top 10 optimized parameter sets (Supplementary file 1a, Figure 2-figure supplement 1) to explore phenomenon of relative sensing in silico. For each parameter set, we simulated the following. The model was first exposed to the background EGF concentration for 50 hr to ensure that all species reached a steady state. The model was subsequently exposed to a step increase in EGF concentration (2-, 3-, 4-, or 6-fold). After the step increase, EGF was kept constant as well. We noted the maximum Akt phosphorylation level at each background concentration and EGF fold-change. For each fit, we obtained a series of maximum pAkt responses across different initial EGF concentrations and for multiple fold changes as well as time integrals of pAkt responses between 0 and 30 min. We combined the predicted relative sensing dose responses at every background EGF level and at every fold by taking the average (and the corresponding standard deviation) across predictions from all 10 best parameter sets. We then plot the resulting dose response as seen in Figure 2c,d of the main text.

Statement of source code availability
All data and source code are available at: https://github.com/dixitpd/FoldChange (Dixit, 2020; copy archived at https://github.com/elifesciences-publications/FoldChange). The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
. Equation A1 describes the dynamics of the concentration [R] of ligand-free receptors. k prod is the rate of delivery of free receptors to cell surface, k 1 is the rate constant of ligand binding to receptors, k -1 is the rate constant of ligand unbinding from the ligand-receptor complex, k 2 is the rate of receptor dimerization, k -2 is the rate of receptor undimerization, k i is the basal internalization rate of non-activated receptors, and finally, k rec is the rate of recycling of internalized unphosphorylated receptors to cell surface. From Equations A1-A8 and the limits imposed by inequalities A14, we calculate the receptor level at steady state with continuous exposure to the ligand by setting the rates of change of individual concentrations to zero. Based on the available knowledge about the biology of EGF receptors and parameters of our model (Equations A14), we neglect the internalization and recycling of ligand-bound monomers and unphosphorylated dimers. Consequently, we neglect degradation of ligand-bound receptor monomers and unphosphorylated receptor dimers. As a result, we assume that at steady state, the fraction of monomeric ligand bound receptors in the endosomes is negligible. We have Solving for steady state Equations A15-A21 when ligand concentration is [L] 0 = u 0 *K d1 , we obtain the formula for the total number of receptors on the cell surface, where we have introduced two aggregate parameters that we denote a and b (see Equations A23a and A23b below). The first term in Equation A22 represents the concentration of monomeric surface receptors (ligand-free [R] and ligand-bound monomeric receptors [LR]) at steady state, the second term represents the steady state concentration of dimeric receptors (unphosphorylated dimers [LR 2 ] and phosphorylated receptors [LR * 2 ]). When the ligand concentration is much lower than the equilibrium dissociation constant; u 0 <<1, we can assume that the majority of cell surface receptors are ligand free as indicated in the second approximation in Equation A22.
Remarkably, even though Equations A1-A8 are governed by more than ten rate parameters, the total concentration of receptors on the surface [R] T at steady state depends only on two composite aggregate parameters a and b and the ligand concentration u 0 . The two aggregate parameters are defined as, a ¼ k dp þ k p k dp The parameter a quantifies the ligand-sensitivity of the signaling system, an increase in the value of a leads to a higher signal sensitivity and an increase in receptor dimerization and phosphorylation. The parameter b quantifies the combination of two biases leading to preferential internalization and degradation of phosphorylated receptors. An increase in the value of b leads to a larger fraction of active (phosphorylated) receptors being internalized and degraded. The first term in A23b is the fraction of dimeric receptors that are phosphorylated. The second term k * i /k i is the relative rate at which phosphorylated receptors are trafficked to the endosomes compared to unphosphorylated receptors. Finally, the third term quantifies the ratio of the two following fractions; k * deg /(k * rec + k * deg ) is the fraction of phosphorylated receptors that are transported from endosomes to be degraded and k deg /(k rec + k deg ) is the fraction of unphosphorylated receptors that are transported from endosomes to be degraded.
Using average values of parameters from our dynamical ODE model fits (Materials and methods and Supplementary files 1, 2, 3), we estimate a~16 and b~40. From here onwards, we use these two values of a and b in our analysis. In Appendix 1-figure 2 we plot the steady state cell surface EGFR levels predicted by Equation A22 as a function of the continuous background ligand stimulation u 0 for these values of a and b. Notably, the receptor levels decrease approximately logarithmically as a function of u 0 over nearly two orders of magnitude in ligand concentration between u 0~1 0 -3.5 -10 À2 .
response is within 80% of the maximum (20% deviation) and so on. Vertical green lines indicate the region of 90% accuracy; the phosphorylation response [LR* 2 ] (scaled) varies within 10% of the maximum. Notably, as a increases, corresponding to a larger fraction of receptors that get phosphorylated, the span of background concentrations over which signaling network exhibits relative sensing increases as well. In panel b), we show [LR* 2 ] (scaled) as a function of background ligand concentration at different values of b. As b increases, the range of background concentrations over which relative sensing holds increases as well. The green arrows show the range of background ligand concentrations over which relative sensing holds for the values of a and b evaluated from the model fits to the data. In summary, the analytical model showed that cells adjust the number of receptors on their plasma membranes in response to background ligand exposures through preferential internalization and subsequent degradation of activated receptors. This effectively allows cells to encode the memory of background ligand exposures on the plasma membrane. The model identified two dimensionless aggregate parameters a and b that dictate the range of background ligand concentrations over which the signaling network can sense relative changes in extracellular ligand concentration. In agreement with experimental data (Figure 3 in main text), the model showed that the relative sensing is observed over an order of magnitude in background ligand exposures and is several orders of magnitude below the equilibrium dissociation constant of the ligand with the receptors. Notably, the model showed that relative sensing was robust to variations in a and b.
Section 4. An analytical model for receptors internalization-based relative sensing for a case of monomer-activated receptors signaling Section 4.1 In this section, we show that relative sensing of extracellular ligand concentration can also occur for signaling cascades where receptor phosphorylation (activation) is initiated by ligand-bound receptor monomers instead of dimers. Therefore, receptor dimerization is not a necessary condition for relative sensing. The derivation of analytical expressions for the monomeric case follows the same general logic used in Sections 1 and 2 above.