1. Cell Biology
  2. Computational and Systems Biology
Download icon

Signaling pathways as linear transmitters

  1. Harry Nunns  Is a corresponding author
  2. Lea Goentoro  Is a corresponding author
  1. California Institute of Technology, United States
Research Article
  • Cited 1
  • Views 1,610
  • Annotations
Cite this article as: eLife 2018;7:e33617 doi: 10.7554/eLife.33617

Abstract

One challenge in biology is to make sense of the complexity of biological networks. A good system to approach this is signaling pathways, whose well-characterized molecular details allow us to relate the internal processes of each pathway to their input-output behavior. In this study, we analyzed mathematical models of three metazoan signaling pathways: the canonical Wnt, MAPK/ERK, and Tgfβ pathways. We find an unexpected convergence: the three pathways behave in some physiological contexts as linear signal transmitters. Testing the results experimentally, we present direct measurements of linear input-output behavior in the Wnt and ERK pathways. Analytics from each model further reveal that linearity arises through different means in each pathway, which we tested experimentally in the Wnt and ERK pathways. Linearity is a desired property in engineering where it facilitates fidelity and superposition in signal transmission. Our findings illustrate how cells tune different complex networks to converge on the same behavior.

https://doi.org/10.7554/eLife.33617.001

Introduction

Cells must continually sense, interpret, and respond to their environment. This is orchestrated by signaling pathways: networks of multiple proteins that transmit signals and initiate cellular response. Signaling pathways are critical to animal development and physiology, and yet there are fewer than 20 classes of metazoan signaling pathways (Gerhart, 1999). These signaling pathways evolved prior to the Cambrian and remain highly conserved across animal phyla (Gerhart, 1999; Pires-daSilva and Sommer, 2003). Each signaling pathway, therefore, governs a wide range of cellular events, both within and across organisms.

Insights into the versatility of signaling pathways may be gleaned from pathway architectures. Indeed, distinct architectural features define each pathway. Studies over the past several decades have revealed distinct signaling capabilities that arise from pathway architecture, for example, all-or-none response in the MAPK/ERK pathway (Huang and Ferrell, 1996; Ferrell and Machleder, 1998), oscillations in the NFκB pathway (Hoffmann et al., 2002), or asymmetrical cell signaling in the Notch/Delta pathway (Sprinzak et al., 2010). Alternatively, analysis of pathway architectures may also reveal shared signaling capabilities that emerge from the distinct architectures, pointing to a fundamental property that pathways have converged upon despite their separate evolutionary trajectories. In this study, we sought to identify shared properties between conserved signaling pathways.

To this end, we examined three signaling pathways, the canonical Wnt, ERK and Tgfβ pathways. These pathways are activated by an extracellular ligand binding to a membrane receptor (Figure 1A). The ligand-receptor activation initiates a series of biochemical reactions within the cell, culminating in a buildup of transcriptional regulator, which regulates transcription of broad gene targets. Since the ligand-receptor module is relatively plastic across organisms (e.g. flies have one EGF receptor whereas humans have four [Citri et al., 2003]), we focused on the conserved core pathway (Figure 1A). We define the input to the core pathway as the ligand-receptor activation, and the output as the level of transcriptional regulator.

The Wnt, ERK, and Tgfβ pathways transmit input using different core transmission architecture.

(A) Signaling pathways transmit inputs from ligand-receptor interaction to a change in output, the level of transcriptional regulator (white circle). (B-D) The core pathway for each metazoan signaling pathway is defined by distinct architectural features. In the Wnt pathway (B), the output is regulated by a futile cycle of continual synthesis and rapid degradation. In the ERK pathway (C), the output is regulated by a kinase cascade coupled to negative feedback. In the Tgfβ pathway (D), the output is regulated through continual nucleocytoplasmic shuttling.

https://doi.org/10.7554/eLife.33617.002

The Wnt, ERK, and Tgfβ pathways transmit input using different core transmission architecture (Figure 1B–D). In the Wnt pathway, signal transmission is characterized by a futile cycle of synthesis and rapid degradation (Kimelman and Xu, 2006; Saito-Diaz et al., 2013; Hoppler and Moon, 2014). We use the term futile cycle to highlight that β-catenin is continually synthesized only to be quickly targeted for degradation and kept at low concentration, as opposed to, for instance, being synthesized only as needed. Ligand-receptor input diminishes the degradation arm of this cycle, leading to accumulation of β-catenin output (Kimelman and Xu, 2006; Stamos and Weis, 2013; Nusse and Clevers, 2017). In the ERK pathway, signal transmission is characterized by a cascade of phosphorylation events coupled to feedbacks, leading to an increase in phosphorylated ERK output (Kolch, 2005; Yoon and Seger, 2006; Avraham and Yarden, 2011; Lake et al., 2016). Finally, signal transmission in the Tgfβ pathway is characterized by continual nucleocytoplasmic protein shuttling (Inman et al., 2002; Nicolás et al., 2004; Xu and Massagué, 2004; Schmierer and Hill, 2005; Massagué et al., 2005). Ligand-receptor input effectively increases the rate of nuclear import, leading to an increase in output, the nuclear Smad complex (Schmierer et al., 2008).

Importantly for our approach, the architectures of the three pathways are captured by mathematical models that have been refined by years of experiments. Although by no means complete, the mathematical models have track records of success in predicting systems-level behaviors across multiple biological systems. For instance, the Wnt model (Lee et al., 2003) captures the dynamics of destruction complex well enough as to enable prediction of robustness in fold-change response (Goentoro and Kirschner, 2009) and the differential roles of the two scaffolds in the pathway (Lee et al., 2003); the ERK model (Huang and Ferrell, 1996; Ferrell and Bhatt, 1997; Schoeberl et al., 2002; Sturm et al., 2010) captures the ultrasensitivity in the phosphorylation cascade (Huang and Ferrell, 1996); and the Tgfβ model (Schmierer et al., 2008) reveals the roles of nucleocytoplasmic shuttling in transducing the duration and intensity of ligand stimulation (Schmierer et al., 2008).

We studied these mathematical models to identify what, if any, behaviors converge across pathways. The Wnt (Lee et al., 2003), ERK (Sturm et al., 2010), and Tgfβ (Schmierer et al., 2008) models consist of 7, 26, and 10 coupled, nonlinear ODEs, respectively, with 22, 46, and 13 parameters. Because of their large sizes, they are typically solved numerically to simulate experimental observations and generate new predictions. However, for the questions posed here, we found that numerical simulations are not sufficient. Rather, we needed analytics to uncover exactly how the pathway behaviors depend on the underlying biochemical processes. While we previously derived an analytical solution to the Wnt pathway (Goentoro and Kirschner, 2009), analytical treatment of the Tgfβ and ERK pathways has not been attempted due to the complex, nonlinear equations involved. To address this problem, we employed various analytical techniques, including graph theory-based variable elimination and dimensional analysis, to derive analytical or semi-analytical solutions to the steady-state output of each pathway. Our analysis, along with subsequent experimental verification, reveals a striking convergence across the Wnt, Tgfβ, and ERK pathways: cells operate in the parameter regime where the complex, nonlinear interactions in each pathway give rise to linear signal transmission.

Results

Mathematical analysis identifies the Wnt, ERK, and Tgfβ pathway as linear transmitters

We began our analysis using established models of the Wnt (Lee et al., 2003), ERK (Sturm et al., 2010), and Tgfβ (Schmierer et al., 2008) pathways. These models capture the salient features of each pathway, and include biochemical details such as synthesis, degradation, binding, dissociation and post-translational modifications. In all the models, biochemical parameters have been directly measured or fitted to kinetic measurements from cell, embryo or extract systems. Numerical simulation of each model has predicted a wide range of pathway behaviors over the years (e.g. Wnt refs. [Lee et al., 2003; Goentoro and Kirschner, 2009; Hernández et al., 2012]; ERK refs. [Huang and Ferrell, 1996; Ferrell and Machleder, 1998; Schoeberl et al., 2002; Sturm et al., 2010; Fritsche-Guenther et al., 2011]; Tgfβ refs. [Schmierer et al., 2008; González-Pérez et al., 2011; Andrieux et al., 2012; Vizán et al., 2013; Wang et al., 2014]). Below, we describe our analysis of each pathway and the unifying behavior that emerges from all three pathways.

Canonical Wnt pathway

In this pathway, cells sense ligand-receptor input by monitoring β-catenin protein (Kimelman and Xu, 2006; Stamos and Weis, 2013; Nusse and Clevers, 2017; MacDonald et al., 2009; Clevers and Nusse, 2012). β-catenin is continually synthesized and rapidly degraded by a large destruction complex, comprised of multiple proteins including APC, Axin, and GSK3β. The destruction complex binds and phosphorylates β-catenin, tagging it for degradation by the ubiquitin/proteosome machinery (Kimelman and Xu, 2006; Stamos and Weis, 2013). Wnt ligands, through binding to Frizzled and LRP receptors, inhibit the destruction complex, leading to accumulation of β-catenin. β-catenin then regulates the expression of broad target genes (Stamos and Weis, 2013; Nusse and Clevers, 2017).

The model of the Wnt pathway (Figure 2A) was published in 2003 by a collaboration between the Kirschner and Heinrich labs (Lee et al., 2003). The Wnt model consists of seven nonlinear differential equations and 22 parameters. Applying dimensional analysis, we previously derived the analytical solution to β-catenin concentration at steady-state (Goentoro and Kirschner, 2009):

(1) [βcat]ss=K171γ+αu2(1+4γ(1γ+αu2)1)
(2) α=k4k6k9v14GSK3totAPCtotk5k_6K7K8k13k15
(3) γ=v12k13K17
Figure 2 with 7 supplements see all
The Wnt, ERK, and Tgfβ pathways are linear signal transmitters.

(A-C) Network diagrams of the signaling pathways. The Tgfβ diagram is modified from Schmierer et al. (2008). In the network diagram in A, DC refers to the β-catenin destruction complex. Below the network diagrams: the parameter groups and linearity equations we analytically derived in this study. Parameter groups and input functions are color-coded to the corresponding reactions in the network diagrams. Parameters that do not appear in the parameter groups either drop out due to irreversible reaction steps (such as k10 and k11 in the Wnt pathway) or negligible (as indicated by ellipses). (D-F) Our analysis reveals that in physiologically relevant parameter values, these pathways generate a linear input-output relationship. The outputs are β-catenin, dpERK, and nuclear Smad complex for the Wnt, ERK, and Tgfβ pathway, respectively. The input functions u describe the effect of ligand-receptor interactions on the core pathway. Specifically: u(Wnt) is the rate by which Dishevelled/Dvl inhibits the destruction complex upon Wnt ligand activation, where k3 and k6 are defined in the figure and [Dvl]a is the concentration Wnt-activated Dishevelled (see Equations A15); u(EGF) is concentration of EGF-activated Ras (Ras-GTP); and u(Tgfβ) is the fraction of Tgfβ -activated receptors. Red and blue lines, respectively: analytical and numerical solutions with measured parameters (plotted against the left y-axis). Grey line: examples of numerical solutions outside measured parameters (plotted against the right y-axis).

https://doi.org/10.7554/eLife.33617.003

where the input function u=u(Wnt) is the rate of inhibition of the destruction complex (DC) via Dishevelled/Dvl, a function of ligand-receptor activation. As illustrated in Figure 2A, Ki’s are equilibrium dissociation constants, ki’s are rate constants, and vi’s are synthesis rates. α and γ in Equation 1 are dimensionless parameter groups defined in Equations 2 and 3: α characterizes β-catenin degradation by the destruction complex, and γ characterizes the extent to which β-catenin binds to APC independently of the destruction complex.

Equation 1 demonstrates that, in general, β-catenin concentration is a nonlinear function of the input u. Many parameters of the model were directly measured in Xenopus extracts, and the remaining calculated from measurements in the same system (Appendix 1—table 1). In this study, we examined how the analytical solution (Equation 1) behaves with these measured parameters. The measured parameters (Appendix 1—table 1) indicate that α~66, γ~1.4, and for maximal stimulation, u~6.0. The large α reflects how β-catenin stability is primarily dictated by the destruction complex, that is, α/u1 means that non-Axin-dependent degradation is minimal, and α/uγ means that the positive feedback from sequestration by APC is minimal. Indeed, the rapid action of the destruction complex in the Wnt pathway is a recurring observation across biological systems (Kimelman and Xu, 2006; Saito-Diaz et al., 2013; Hoppler and Moon, 2014). With α/u1+γ, Equation 1 simplifies to

(4) [βcat]ssK17γαu

with detailed derivations presented in Appendix 1. Therefore, within physiologically relevant parameter values, the steady-state β-catenin concentration becomes a linear function of the input u (red line, Figure 2D). The linear input-output relationship holds for the entire dynamic range of the model, until the system saturates at maximal stimulation (u6.0). We confirmed that the numerical solution of the full model matches the analytical solution in Equation 4 (blue line, Figure 2D), and that the response becomes nonlinear when α is decreased, breaking the requirement α/u1+γ (grey line, Figure 2D).

Source codes for the numerical simulations in Figure 2D–F (grey and black lines) are available in Figure 2—source code 1.

ERK pathway

The unexpected linearity that emerges from the model of the Wnt pathway prompted us to wonder if such simplicity may be found in other pathways. Strikingly, we observed the same linearity in the ERK and Tgfβ pathways. In the ERK pathway (Figure 2B), ligand-receptor input is transmitted via a cascade of protein phosphorylation (Kolch, 2005; Yoon and Seger, 2006). In particular, ligand-receptor interactions activate Ras, which leads to membrane recruitment and phosphorylation of Raf. Phosphorylated Raf subsequently doubly phosphorylates MEK, which in turn doubly phosphorylates ERK (Kolch, 2005). Doubly-phosphorylated ERK (dpERK) is a transcriptional regulator that affects a broad array of genes (Yoon and Seger, 2006). The multi-step topology of the kinase cascade, combined with distributive phosphorylation of each kinase, gives rise to ultrasensitivity – first demonstrated in the seminal work by the Ferrell lab (Huang and Ferrell, 1996; Ferrell and Machleder, 1998). In other contexts, the pathway also exhibits a graded response (Whitehurst et al., 2004; Mackeigan et al., 2005; Cohen-Saidon et al., 2009; Ahmed et al., 2014) that is thought to arise from the incorporation of negative feedbacks (Lake et al., 2016), one of which is the inhibition of Raf by dpERK through hyper-phosphorylation of serine residues (Sturm et al., 2010; Dougherty et al., 2005; Hekman et al., 2005).

The ERK model (Sturm et al., 2010) is the product of more than two decades of refinement (Huang and Ferrell, 1996; Ferrell and Machleder, 1998; Schoeberl et al., 2002; Sturm et al., 2010; Fritsche-Guenther et al., 2011). The model, which captures ultrasensitivity and Raf feedback, consists of 26 differential equations and 46 parameters. To derive an analytical expression for the ERK pathway, we used a variable elimination technique developed for networks of mass action kinetics (Feliu and Wiuf, 2012). The technique utilizes an algebraic framework, linear elimination of variables, and mass conservation laws to parameterize steady-state in terms of core variables (described in Appendix 1). We derived an analytical relationship between the steady-state output of the pathway [dpERK]ss and the input to the phosphorylation cascade u:

(5) [dpERK]ss=αβ(Raftot[pRaf]ss)1γαuδβ
(6) α=k3(k8+kb7)k7[P1]ssk8+
(7) β=k25(k30+kb29+k29[P4]ss)k29[P4]ssk30+
(8) γ=k3(k8+kb7)k9[MEK]ssk7[P1]ssk8k10+
(9) δ=k26+kb25k26+

Detailed derivations of Equation 5 are presented in Appendix 1. The input u=u(EGF) in Equation 5 is the concentration of active Ras, which is activated via GTP loading at the ligand-receptor complex (Kolch, 2005). The parameter groups α, β, γ, and δ in Equation 5 are defined in Equations 6–9, where the ellipses indicate additional small terms (expanded in Appendix 1). The relative magnitudes of α, β, γ, and δ indicate how the Raf pool partitions during signaling (Equations A21, A29A31). The dimensionless group αu relates to the amount of free, phosphorylated Raf (α, blue-shaded in Figure 2B), β[dpERK]ss describes the amount of Raf inhibited through negative feedback by dpERK (β, red-shaded in Figure 2B), δ relates to the amount of unphosphorylated (δ, blue-shaded in Figure 2B), and γu relates to the amount of phosphorylated Raf bound to other proteins (e.g. to MEK, brown-shaded in Figure 2B). Equation 5 is not a closed solution, as it includes the term pRafss, and there are variables included in parameter groups α, β, γ. We confirmed that the parameter groups remain constant over the course of signaling (within 10%, Figure 2—figure supplement 1), justifying treating the latter variables as parameters.

Next, we considered how the analytical expression (Equation 5) behaves within a specific parameter regime observed in experiments. First, experiments in several mammalian cell systems have shown that feedback is strong, such that a significant fraction of the Raf pool is inhibited (Fritsche-Guenther et al., 2011; Dougherty et al., 2005). This means that βdpERKss~α+γu+δ. Second, as has been observed in multiple contexts ([Huang and Ferrell, 1996; Ferrell and Machleder, 1998; Schoeberl et al., 2002; Sturm et al., 2010] Appendix 1—table 2), ERK phosphorylation is ultrasensitive to the amount of pRaf (the ultrasensitive cascade is shaded green in Figure 2B). Denoting K as the relative change of dpERKss with respect to pRafss, ultrasensitivity entails that K1. In this range, small changes in pRaf level have very large effects on dpERK level (e.g., in model simulations, a 30% change in pRaf level results in a 900% change in dpERK level, Figure 2—figure supplement 1). We find analytically that in the parameter regime where βdpERKss~α+γu+δ and K1, the negative feedback holds the level of pRaf constant (pRafssRs, details in Appendix 1). With these two features, strong negative feedback and ultrasensitivity, dpERK becomes a linear function of the input u:

(10) [dpERK]ssαβRaftotRsuδβ

The full derivation is given in Appendix 1, and includes a toy model to illustrate the intuition for how ultrasensitivity combines with negative feedback to produce linearity. Equation 10 is plotted in Figure 2E (red line). We confirmed that the numerical solution of the full model matches the analytics in Equation 10, and becomes nonlinear when the negative feedback is weakened (grey line, Figure 2E). Although the analytical expression describes up until 50% of ERK activation, we verified numerically that the predicted linearity extends to 93% of ERK activation (Figure 2—figure supplement 2).

The linearity derived here applies across different dynamic ERK responses. The model we analyzed gives a sustained dpERK response. In some contexts, however, the ERK pathway shows a pulsatile response, which has been attributed to receptor desensitization (Schoeberl et al., 2002). Using a larger model that includes details of receptor desensitization (Schoeberl et al., 2002), we numerically verified that the linearity holds for pulsatile responses - that is, the peak level of dpERK increases linearly with the peak level of u (Figure 2—figure supplement 1).

Tgfβ pathway

Finally, we examined signal transduction within the Tgfβ pathway (Figure 2C). In the Tgfβ pathway, input from ligand-receptor interactions is transmitted by the Smad proteins. There are several classes of Smad proteins, including the receptor-regulated Smads (R-Smads) and the common Smad (co-Smad or Smad4) (Massagué et al., 2005). Ligand-activated receptors phosphorylate R-Smads. Phosphorylated R-Smads bind to the co-Smad, and shuttle into the nucleus and regulate broad target genes. In the nucleus, the Smad complex dissociates and R-Smads are constitutively de-phosphorylated and shuttled out to the cytoplasm, where the cycle of phosphorylation and complex formation begins again (Schmierer et al., 2008). This dynamic translocation in and out of the nucleus forms a continual nucleocytoplasmic shuttling of Smads, a known integral feature of the Tgfβ pathway (Inman et al., 2002; Nicolás et al., 2004; Xu and Massagué, 2004; Schmierer and Hill, 2005).

The Tgfβ model (Schmierer et al., 2008) was published in 2008 by the Hill lab, and consists of 10 differential equations and 13 parameters. Even though the model was fitted to R-Smad2 data, the general architecture of signal transmission is conserved across all five R-Smads (Massagué et al., 2005; Schmierer et al., 2008). Using the variable elimination technique described before (Feliu and Wiuf, 2012), we derived an analytical expression of the steady-state concentration of Smad complex in the nucleus:

(11) [S24n]ss=aαu(α+γ)u+βS2tot
(12) α=a(kon[S4n]ss+akex2)koff+
(13) β=PPasekdephoskphosRtotkex2akex2+kin2+
(14) γ=a(akex2+PPasekdephos)(1akex2+1CIFkin2)+

In Equation 11, the input function u=u(Tgfβ) is the active fraction of Tgfβ receptors. The parameter a is the nucleocytoplasmic volume ratio. The dimensionless parameter groups α, β, and γ in Equation 11 are defined in Equations 12–14, where the ellipses indicate additional small terms (expanded in Appendix 1). α, β, and γ describe how the Smad2 pool partitions during signaling (Equations A44, A50, A51): αu relates to the amount of nuclear Smad complex (α, blue-shaded in Figure 2C, captures the parameters related to complex formation and translocation to the nucleus), β relates to the amount of free, unphosphorylated Smad2 (β, red-shaded in Figure 2C, captures the parameters related to complex dissociation and translocation to the cytoplasm), and γu loosely relates to the remaining Smad2 pool (γ is brown-shaded in Figure 2C). Phosphorylated Smad2 quickly forms complex (Lagna et al., 1996), so β essentially corresponds to total monomeric Smad2. Finally, Equation 11 is not a closed solution, since variable S4nss appears in α. We numerically tested that it is constant within 2% for non-saturating inputs (Figure 2—figure supplement 3), justifying treating it as a parameter.

As in the Wnt and ERK pathway, the analytical expression for nuclear Smad complex (Equation 11) allows us to see that the behavior dramatically simplifies with parameters observed in experiment. We consider the case for non-saturating inputs (u0.1). Protein concentrations in the Tgfβ model were measured in human keratinocyte cells and the rate constants fitted to kinetic data measured in the cells (Schmierer et al., 2008). With the measured parameters (Appendix 1—table 3), we find that β~46, αu~1.5, and γu~0.7. In this parameter regime, once Smad2 is imported to the nucleus, it is rapidly dephosphorylated and exported. Dynamic Smad2 translocation maintains monomeric Smad2 in excess to Smad complex (βα+γu). and forms the continual nucleocytoplasmic shuttling that is characteristic of the Tgfβ pathway. Even under maximal Tgfβ stimulation, it has been estimated that phosphorylated Smad2 comprises only 36% of the Smad2 pool (Schmierer and Hill, 2005; Gao et al., 2009). With βα+γu, the first term in the denominator of Equation 11 is small, and concentration of nuclear Smad complex becomes a linear function of input:

(15) [S24n]ssaαS2totβu

Equation 15 is plotted in Figure 2F (red line), and we confirmed that numerical simulations recapitulates Equation 15 (blue line, Figure 2F). Although the analytical solution is valid only for small values of u, we numerically verified that the predicted linearity holds for the entire range of input u (from 0 to 1, Figure 2—figure supplement 2). We confirmed that the pathway becomes nonlinear when the R-Smad phosphatase is inhibited such that β~α+γu (grey line, Figure 2F). While the model analyzed here gives a sustained Smad response, we verified numerically that the linearity holds for a larger model that includes receptor desensitization and gives a pulsatile Smad response (Figure 2—figure supplement 3) (Vizán et al., 2013).

Linearity in the Wnt and ERK pathways was observed experimentally

Analytical expressions for the Wnt, ERK, and Tgfβ pathways reveal that the three pathways behave as linear signal transmitters within parameter regimes measured in cells. To confirm the linearity, we directly measured the input-output relationships in human cell lines. We focused our efforts on the Wnt and ERK pathways, since we are limited by available antibodies in the Tgfβ pathway.

To analyze the canonical Wnt pathway, we performed quantitative Western blot measurements in RKO cells, a model system for Wnt signaling. To track the input, we measured the level of phosphorylated LRP5/6 receptors (on Ser1490), which increases within minutes of ligand-receptor complex formation (Tamai et al., 2004). To track the output, we measured the level of β-catenin. We confirmed that the level of phosphorylated LRP5/6 and β-catenin increase upon Wnt simulation and reach steady-state within 6 hr (Figure 3—figure supplement 1). Accordingly, all subsequent measurements were done at 6 hr after Wnt stimulation.

To measure the input-output relationship in the Wnt pathway, we treated RKO cells with varying doses of purified Wnt3A and measured how β-catenin (output) correlates with phosphorylated LRP (input). As shown in Figure 3A, the level of β-catenin increases linearly with the level of phosphorylated LRP. The linearity persists until saturation of the input, defined as 90% of maximal phosphorylated LRP response (blue circles, Figure 3A; Figure 3—figure supplement 2). Notably, at high doses of Wnt3A, β-catenin continues to show incremental activation, despite saturation in phosphorylation of LRP (grey circles, Figure 3A). This can be explained within some findings that, while Frizzled/LRP complex is the primary receptor input in β-catenin activation, β-catenin can be activated independently of LRP (e.g. Rotherham and El Haj, 2015).

Figure 3 with 9 supplements see all
Linearity was observed experimentally in the Wnt and ERK pathways.

(A) Measurements of the input-output relationship in the Wnt pathway. In these experiments, RKO cells were stimulated with 0–1280 ng/mL purified Wnt3A ligand, harvested at 6 hr after ligand stimulation, and lysed for Western blot analyses. Shown on top is a representative Western blot. The data plotted come from seven independent experiments (total N = 66). Each circle indicate the mean intensities of the phospho-LRP5/6 (x-axis) and β-catenin (y-axis) bands for all Western blot biological replicates, and error bars indicate the standard error of the mean. For each gel, we normalize the unstimulated sample (i.e. 0 ng/mL of Wnt3A) to one, and scale the magnitude of the dose response to the average of all gels (described in Materials and methods). The grey line is a least squares regression line, and ρ is the Pearson’s coefficient, where ρ = 1 is a perfect positive linear correlation. (B) Measurements of the input-output relationship in the ERK pathway. In these experiments, H1299 cells were stimulated with 0–50 ng/mL purified EGF ligand, harvested at 5 min after ligand stimulation, and lysed for Western blot analyses. Shown on top is a representative Western blot. The data plotted here come from five independent experiments (total N = 30). Each circle indicates the mean intensities of dpERK1/2 bands across Western blot biological replicates, and the error bars indicate standard error of the mean. Single replicates are plotted without error bars. All data is plotted relative to unstimulated sample. The grey line is a least squares regression line, and r2 is the coefficient of correlation where r2 = 1 is a perfect linear correlation. (C) As in (A), except that cells were treated with 1 μM CHIR99021 (detailed in Materials and methods). The data plotted here come from five independent experiments (total N = 59). The grey line is a least squares regression, and ρ is the Pearson’s coefficient, where ρ = 1 is a perfect positive linear correlation. Shown in the subplot are the same least squares regression line (solid line), overlaid with the model prediction (dashed line). (D) As in (B), but measurements were performed in H1299 cells expressing mutant Raf S289/296/301A. The data plotted here come from three independent experiments (total N = 15). The grey line is a fit using the ERK model. We first fitted the gain of the model to the data (i.e. the y-range), and afterward, varied the strength of dpERK feedback (k25) to find the best fit. We used the weighted Akaike Information Criterion, w(AICc), to verify that the nonlinear fit from the ERK model outperforms a linear least squares fit (see Materials and methods). 0 < w(AICc) < 1, with higher w(AICc) indicates better performance by the non-linear fit. In all figures, linearity was additionally assessed using the least absolute deviations, L1-norm (see Methods). L1-norm can range from 0 to 0.5, with L1-norm < 0.1 indicate a linear relationship. Blue vs grey circles in each figure are explained in the main text. Source files of all Western blot gel images and numerical quantitation data are available in Figure 3—source data 1.

https://doi.org/10.7554/eLife.33617.012

Consistent with the mathematical analysis, we observed in RKO cells that the Wnt pathway behaves as a linear transmitter throughout the dynamic range of the input. As a control that is expected from the Michaelis-Menten kinetics that describe ligand binding in the model, we confirmed that the linearity does not extend upstream to Wnt dose: both phospho-LRP5/6 and β-catenin show nonlinear response to Wnt dose (Figure 3—figure supplement 2). Therefore, in the Wnt pathway, a nonlinear ligand-receptor processing step is followed by linear signal transmission through the core intracellular pathway.

Next, to measure the input-output relationship in the ERK pathway, we performed quantitative Western blots in H1299 cells, one of the model systems used in the field. Linearity in the ERK pathway has been suggested in different parts of the pathway, e.g. Knauer et al. (1984) used experimental and modeling analyses to infer linearity between receptor occupancy and the downstream cellular proliferation; Oyarzún et al. (2014) suggests linearity in ligand-receptor processing. Here, we specifically probe linearity in the core transmission step of the pathway. Detecting the input level, EGF-activated Ras GTP, requires a pull down step that makes it less quantifiable. Therefore, motivated by Oyarzún et al. (2014), we tested EGF ligand itself as the input. To track the output, we measured the level of doubly-phosphorylated ERK1/2 (on Thr202/Tyr204), dpERK. We first characterized the kinetics of response: dpERK peaks 5 min after EGF stimulation (Figure 3—figure supplement 3), and saturates at 4 ng/ml EGF (grey circles, Figure 3B). Accordingly, all subsequent measurements were performed at 5 min after EGF stimulation, and linearity was assessed over the input range of 0–4 ng/mL EGF (blue circles, Figure 3B).

We observed linearity in the input-output relationship of the ERK pathway, with the level of dpERK increasing linearly with EGF dose (Figure 3B). The linearity holds throughout the dynamic range of the system, over at least 12-fold activation of dpERK. As the ERK pathway is sometime observed to show bimodal response that would be masked by bulk measurements, we confirmed that the H1299 cells indeed show to graded dpERK response in single-cell level (Figure 3—figure supplement 4), in agreement with a previous single-cell, live imaging study (Cohen-Saidon et al., 2009). Therefore, as in the Wnt pathway, signals are transmitted linearly in the ERK pathway throughout the dynamic range of the cell. Moreover, the linearity in the ERK pathway is more extensive than in the Wnt pathway, as linearity extends all the way upstream, such that the level of dpERK directly reflects the dose of extracellular EGF ligand.

Linearity in the Wnt and ERK pathways is modulated by perturbation to parameters

Finally, the analytical expressions we derived in this study not only reveal linear signal transmission, but also the mechanisms by which it arises. In the model of the Wnt pathway, linear transmission occurs due to the futile cycle of β-catenin, in the parameter regime where β-catenin is continually synthesized and rapidly degraded (i.e. α/u1+γ). This regime is not infinite: for instance, a ten-fold decrease in α (e.g. by inhibiting the destruction complex) will break the futile cycle (grey line, Figure 2D).

To test if the futile cycle is indeed required for linear signal transmission, we inhibited the destruction complex using CHIR99021, an inhibitor of GSK3β kinase. As before, we measured the input-output relationship, β-catenin vs. phospho-LRP5/6 level, up to 90% of maximal phospho-LRP5/6 input (blue circles, Figure 3C). As expected, we found that inhibiting the destruction complex (decreasing α in the model) reduced the range of linearity. The non-treated cells (blue circles, Figure 3A) exhibit a linear input-output relationship over a 4.4-fold range of LRP input, whereas the CHIR-treated cells show a linear input-output relationship over only a 2.8-fold range of LRP input (blue circles, Figure 3C).

Further, our measurements also reveal an unexpected feature of the Wnt pathway. In the model, inhibiting GSK3β causes β-catenin response to become nonlinear for larger inputs (dashed line, Figure 3C subplot). In CHIR-treated RKO cells, however, this nonlinearity cannot be reached, as the maximal amount of phosphorylated LRP (input) is reduced by 50% (grey circles, Figure 3; Figure 3—figure supplement 2), consistent with the dual-function of GSK3β identified by Zeng et al. (2005); Zeng et al. (2008) in phosphorylating β-catenin for degradation as well as phosphorylation LRP for activation. Incorporating this dual-role of GSK3β into the model, we found that this expanded model can indeed recapitulate the data (Figure 2—figure supplement 4). Therefore, our data indicate two findings: first, that inhibiting GSK3β reduces the range of linear input-output behavior in the Wnt pathway, as predicted by our analytics, and second, that GSK3β co-regulation of β-catenin and LRP unexpectedly constrains the system within the linear regime.

Next, we examine the requirements for linearity in the ERK pathway. Equation 10 reveals that linearity in the ERK pathway depends upon the coupling of strong nonlinearities – ultrasensitivity and negative feedback. As in the Wnt pathway, this regime is not infinite, for example, decreasing the strength of feedback β enables the system to exit the ultrasensitive regime, and therefore reduces linearity (grey line, Figure 2E).

To test this requirement, we examined the effects of weakening the negative feedback. We created a stable H1299 cell line expressing Raf S289/296/301A, a Raf-1 mutant in which three serine residues that are phosphorylated by dpERK are mutated to alanine (Dougherty et al., 2005; Hekman et al., 2005). Assessing the dynamic range of the input as before (0–4 ng/mL EGF), we now found that dpERK responds nonlinearly to EGF dose (blue circles, Figure 3D), consistent with model predictions (grey line, Figure 3D). As a control, we found that overexpressing WT Raf-1 to a similar level does not perturb linearity (experiments, Figure 3—figure supplement 5; modeling, Figure 2—figure supplement 1). Lastly, mutating all five direct ERK feedback sites on Raf-1 to alanine had a similar effect to Raf S289/296/301A (Figure 3—figure supplement 6). Our results support the model requirement that strong negative feedback is critical to linear signal transmission in the ERK pathway.

Discussion

Our study suggests that the canonical Wnt pathway, the ERK pathway, and the Tgfβ pathway have converged upon a shared strategy of linear signal transmission. Our mathematical analysis reveals that, despite their distinct architectures, the three signaling pathways behave in some physiological contexts as linear transmitters. Not only is linearity is predicted within measured parameter regimes, the analysis shows that linearity is a property of the systems that occurs through a considerable range of parameters (Figure 2—figure supplements 5 and 6). We then showed direct measurements of the linear input-output relationship in the canonical Wnt and ERK pathway.

It would be interesting to further probe the generality of linear signal transmission. Linear behavior requires that single cells responds to ligand in a graded manner. Although there are reports of oscillatory or bimodality in signaling pathways, there are also multiple observations across biological contexts of single cells responding to ligand in a graded manner (Appendix 1—table 4). Besides the systems analyzed here, NF-κB is another signaling pathway that has been modeled rigorously (Hoffmann et al., 2002; Ashall et al., 2009; Lee et al., 2014). Numerical simulations of a well-established NF-κB model (Ashall et al., 2009) over the range of nuclear NF-κB translocation observed in human epithelial cells (Lee et al., 2014) reveal that the peak of the nuclear NF-κB pulse correlates linearly with ligand concentration (Figure 2—figure supplement 7). Finally, linearity extends beyond metazoan signaling pathways. In the yeast pheromone sensing pathway, a homolog of the ERK cascade, transcriptional output correlates linearly with receptor occupancy (Yu et al., 2008). The linearity is mediated by negative feedback by Fus3 acting on Sst2, a feedback that is not conserved in the mammalian ERK system. These further argue for linear signal transmission as a convergent property across independently evolving signaling pathways, as well as between conserved pathways that diverged 1.5 billion years ago.

What are potential advantages to linear signal transmission? Linearity is a feature of many engineering systems, where it serves several practical purposes. In particular, linear signal transmission enables the superposition of multiple signals, where the output of two simultaneous inputs is equal to the sum of the outputs for each input separately. Superposition enables multiple, dynamic signals to be faithfully transmitted and processed independently. Thus, for instance, linearity enables people to listen to a phone call and interpret speech amongst background noise, and allows a car radio to tune into one station out of multiple broadcasting on separate carrier frequencies. Notably, linearity is also a desired goal in synthetic biology, where it is often implemented using negative feedback (Nevozhay et al., 2009; Del Vecchio et al., 2016). Analogous to engineered circuits, linearity in biological signaling pathways may facilitate multiplexing inputs into a single pathway (Figure 4A).

Benefits of linearity.

(A) Linearity enables multiplexing of inputs to a signaling pathway. Multiplexed signals can be independently decoded downstream, and therefore regulate distinct transcriptional events. (B) Illustration for how linearity between the receptor occupancy and downstream outputs gives rise to dose-response alignment (Andrews et al., 2016). (C) Linearity can produce fold-changes in output that are robust to variation in cellular parameters. To illustrate this, we added lognormal noise (0.1 CV) to all parameters of the Wnt model, and simulated the level of β-catenin before and after Wnt stimulation (blue circles). As long as the model operates in the regime of linear signal transmission (i.e. Y=au, where Y is output, u is input, and a is a scalar that is a function of parameters), variation in parameters affects stimulated and basal level of β-catenin equally, and we get a constant fold change in β-catenin (i.e. red line, where FC=Ystimulated/Ybasal is independent of parameter variations).

https://doi.org/10.7554/eLife.33617.023

A second benefit is that linearity might underlie two phenomena that are increasingly found across signaling pathways. First, a linear transmitter naturally gives rise to dose-response alignment (Andrews et al., 2016), where one or more downstream responses of a pathway closely follows the fraction of occupied receptor (Figure 4B). Dose response alignment appears in many biological systems and is thought to improve the fidelity of information transfer through signaling pathways (Oyarzún et al., 2014; Yu et al., 2008; Andrews et al., 2016; Becker et al., 2010). Second, linearity facilitates fold change detection, where cells sense fold changes in signal, rather than absolute level, to buffer cellular noise (Goentoro and Kirschner, 2009; Cohen-Saidon et al., 2009; Lee et al., 2014; Thurley et al., 2014; Frick et al., 2017). In linear input-output systems, the stimulated output correlates linearly to the basal output; thus, the fold-change in output is robust to variations in cellular parameters (Figure 4C). Indeed, for the signaling pathways studied here, it has been shown experimentally that the robust outcome of ligand stimulation is the fold-change in the level of transcriptional regulator (Goentoro and Kirschner, 2009; Cohen-Saidon et al., 2009; Lee et al., 2014; Frick et al., 2017). Therefore, selecting for linearity may naturally confer the benefits of superposition, dose-response alignment, and a robust fold-change in output.

Interestingly, unlike synthetic circuits whose linearity is often designed to extend across multiple orders of magnitude (Nevozhay et al., 2009; Nevozhay et al., 2013), the linearity we observed in the three natural pathways extends only one order of magnitude, which is also the dynamic range of the pathways. However, we know that natural pathways can convey inputs varying across multiple orders of magnitude, for example, vision. Thus, an advantage of linearity in natural pathways may be that, in conjunction with fold-change detection at the receptor-level (Olsman and Goentoro, 2016), the system as a whole can continually adapt to a given input, hence maintaining sensitivity to future signals.

Why evolve complexity in signaling pathways only to produce seemingly simple behavior? We offer two thoughts. First, complexity of each pathway might afford tunability, in the sense that parameters can be tuned to produce different behaviors in different contexts. For instance, the ERK pathway produces digital, all-or-none response in some contexts (Huang and Ferrell, 1996), and analog response in others (Whitehurst et al., 2004; Mackeigan et al., 2005). Second - to take an example from engineering - in order to utilize physical processes that are not naturally linear, engineers must implement complex design features to approximate linearity. Similarly, many biochemical processes are inherently nonlinear, meaning that linearity does not arise from a reduction in complexity. Indeed, in each pathway we analyzed here, linearity emerges from complex interactions: a futile cycle in the Wnt pathway, ultrasensitivity coupled to feedback in the ERK pathway, and continual nucleocytoplasmic shuttling in the Tgfβ pathway. Therefore, analogous to engineered systems, complexity in the biochemical pathways we analyzed here might have evolved in part to produce linearity.

Materials and methods

Expression constructs

pBABEpuro-CRAF that contains the wt human Raf-1 clone was a gift from Matthew Meyerson (Addgene plasmid # 51124). Mutant Raf (S289/296/301A) and (S29/289/296/301/642A) were generated using the Q5 site-directed mutagenesis kit (New England Biolabs, E0554S). The mutant and wt Raf-1 were then placed downstream of a CMV promoter.

Cell lines and cell culture

RKO cells (ATCC, CRL-2577) and H1299 cells (ATCC, CRL-5803) were authenticated by STR profiling and supplied by ATCC. RKO cells were cultured at 37°C and 5% (vol/vol) CO2 in DMEM (ThermoFisher Scientific; 11995) supplemented with 10% (vol/vol) FBS (Invitrogen; A13622DJ), 100 U/mL penicillin, 100 μg/mL streptomycin, 0.25 μg/mL amphotericin, and 2 mML-glutamine (Invitrogen). H1299 cells were cultured at 37C and 5% (vol/vol) CO2 in RPMI (ThermoFisher Scientific; 11875) supplemented with 10% (vol/vol) FBS (Invitrogen; A13622DJ), 100 U/mL penicillin, 100 μg/mL streptomycin, 0.25 μg/mL amphotericin, and 2 mML-glutamine (Invitrogen). Both cell lines tested negative for mycoplasma contamination.

Transfection of Raf-1 constructs

H1299 cells were transfected with the mutant and wt Raf-1 constructs using Lipofectamine 3000 (ThermoFisher Scientific, L3000). Stable expression was selected using puromycin at a concentration of 1.5 μg/mL for 2 weeks.

Reagents and antibodies

The following antibodies were purchased from Cell Signaling Technologies: anti-Phospho-p44/42 MAPK (Erk1/2) (Thr202/Tyr204) (E10) Mouse mAb #9106, anti-histone H3 (D1H2) XP Rabbit mAb #4499, anti-c-Raf Antibody #9422, anti-phospho-LRP6 (Ser1490) Antibody #2568, anti-GAPDH (D4C6R) Mouse mAb #97166. Anti-Beta-catenin mouse mAb was purchased from BD Transduction Laboratories (#610153) and anti-GAPDH rabbit antibody was purchased from Abcam (ab9485). The following fluorescent secondary antibodies were purchased from Fisher Scientific: IRDye 800CW Goat anti-Mouse IgG (926–32210) and IRDye 680LT Goat anti-Rabbit IgG (926-68021).

Recombinant human Wnt3A was purchased from Fisher Scientific (5036WN), and recombinant human EGF was purchased from Sigma (E9644). CHIR99021 was purchased from Sigma (SML1046). Halt Protease and Phosphatase Inhibitor Cocktail (100X) was purchased from Fisher Scientific (78440).

CHIR99021 treatment

RKO cells were pre-treated with 1 μM CHIR99021 for 24 hr before adding replacement media containing 1 μM CHIR99021 and Wnt3A for 6 hr.

Cell lysis

RKO cells at 70% confluency were scraped in PBS, pelleted, and snap-frozen, and then thawed in NP-40 lysis buffer containing Halt inhibitor cocktail. Samples were spun down, and the supernatants were transferred to Laemmli sample buffer and boiled. The samples were then run onto a Bolt 4–12% Bis-Tris Plus Gel (Thermofisher, NW04120BOX). H1299 cells at 70% confluence were scraped in NP-40 lysis buffer containing Halt inhibitor cocktail, and further lysed in Laemmli sample buffer. Samples were spun down, and the supernatants were boiled. The samples were then run onto a Novex 4–20% Tris-Glycine Mini Gel (ThermoFisher, XP04200BOX).

Quantitative Western blots

Proteins were transferred onto nitrocellulose membranes, blocked for one hour at room temperature (RT) with blocking buffer (Odyssey Blocking Buffer (TBS) (927–50000) or 5% milk powder in TBS) and stained overnight at 4°C with primary antibody diluted in blocking buffer. The membranes were then stained with fluorescent IR secondary antibodies diluted in blocking buffer for one hour at RT. The fluorescent signal was then imaged using the LiCOR Odyssey Imager and quantified using Odyssey Application software version 3.0. The background-subtracted intensity of the protein bands were normalized to the loading control, GAPDH and/or Histone H3 (for RKO) or Histone H3 (for H1299). These values were then normalized to the reference lanes within each gel, to allow comparison across gels. For β-catenin and phospho-LRP5/6, variation in the fold-activation from experiment to experiment could artificially stretch the data along the x- and y-axis, and introduce artifacts into the relationship between phospho-LRP5/6 and β-catenin. Therefore, for Wnt3A dose responses, the data from each gel was scaled such that the mean of 80 ng/mL and 160 ng/mL samples was equal to the mean across all gels. Finally, for each antibody used in the study, we did careful characterization of the linear range, and verified that our measurement conditions were within the linear range of the antibody. Technical variability of Western blot quantitation. To confirm the effects reported, we verified that quantitation of the same sample loaded in multiple lanes in a gel gives CV < 10%, and quantitation of the same sample across multiple independent gels gives CV < 10% (Figure 3—figure supplement 7). As further control, we verified that normalization with loading control did not produce artificial distortion of the input-output relationship: linearity was observed without normalization in cases where loading was already uniform (Figure 3—figure supplement 8).

L-1 and L2-norm analysis

L1-norm analysis was performed as described in Nevozhay et al. (2013). Briefly, the data is fitted with a cubic Hermite polynomial, and rescaled along the x and y axis to [0, 1]. The L1-norm is computed as the area between the polynomial fit and the diagonal. Linearity is defined in this context as L1-norm < 0.1. L2-norm analysis for Wnt pathway data was performed using a Pearson’s coefficient, and L2-norm analysis for ERK pathway data was performed using the coefficient of correlation, r2.

Akaike information criterion

To score the validity of nonlinear model fits for Figure 3D, we used the bias-corrected Akaike Information Criterion as described in ref. (Spiess and Neumeyer, 2010), which assesses goodness-of-fit and model parsimony. The weighted Aikaike w(AIC) provides a comparison of all considered models, which in our case is the nonlinear ERK pathway model fit and a linear fit, with the higher score indicating a more valid model.

Appendix 1

1. Variable elimination

We use a variable elimination technique from Feliu and Wiuf (2012) to derive analytic expressions for the steady-states of the Tgfβ and ERK pathways. This technique was developed to handle the complexity of large chemical reaction networks. By eliminating variables from the steady-state solution, we can express the steady-state of the system in terms of a smaller subset of variables. This is a useful tool for analyzing the Tgfβ and ERK models, as the steady-state solution consists of a large set of variables, each with a polynomial equation describing its steady-state.

The technique works as follows: if we can identify a cut set within the reaction network, we can reduce the system to a set of first-order homogeneous equations with respect to that cut. This set of equations can then be solved using linear algebra.

A cut is a set of species such that for every reaction involving those species, there is exactly one reactant and one product that falls within that cut. For example, let us consider a network of four interacting species, A, B, C, and D.

Appendix 1—scheme 1
Network of four proteins.
https://doi.org/10.7554/eLife.33617.026

In this network, there is a cut {A,B,C} that contains exactly one product and one reactant for each reaction. We have highlighted this cut in the reaction set:

Appendix 1—scheme 2
Reaction set corresponding to protein network.
https://doi.org/10.7554/eLife.33617.027

The species D cannot belong in the cut, since it appears twice as a reactant in the first reaction. The behavior of this network is described by four differential equations,

(A1) [A].=k1[A][D]2+k3[C]=0
(A2) [B].=k1[A][D]2k2[B]=0
(A3) [C].=k2[B]k3[C]=0
(A4) [D].=2k1[A]D2+k2[B]+k3[C]=0

which are set to zero at steady-state, and two additional conservation equations:

(A5) T1=A+B+[C]
(A6) T2=2B+C+[D]

The variable elimination technique allows us to reduce the steady-state system of equations by four (three equations for the cut set, and one conservation equation). We do this by expressing each member of the cut set as a dependent variable of D, shown below. We utilize the fact that the differential equations for A, B, and C are first-order and homogenous with respect to our cut, and rewrite them in matrix form. We use the subscript ‘ss’ to denote steady-state:

(A7) -k1[D]ss20k3k1[D]ss2-k200k2-k3[A]ss[B]ss[C]ss=0

Feliu and Wiuf (2012) provides a proof of why a cut set guarantees that we can rewrite the corresponding equations in matrix form. It can be understood intuitively from the fact that a cut contains exactly one reactant of each reaction, and therefore each rate is first-order with respect to the cut. Homogeneity also follows from this, since there are no rate terms that do not include members of the cut.

For a complex model, there is no guarantee that we can derive closed-form analytical solutions for steady-state. The matrix formulation and variable elimination technique immediately provides us with a set of solvable variables. The solution to the matrix equation above is:

(A8) [A]ss=ck2k3
(A9) [B]ss=ck1k3Dss2
(A10) [C]ss=ck1k2Dss2

c is a scaling factor not constrained by the matrix equation. With the use of the conservation Equation 5, we can calculate c and express the steady state of all three species solely in terms of the parameters of the network, and [D]ss. For instance, the solution for [C]ss is below.

(A11) [C]ss=k1k2[D]ss2k2k3+k1[D]ss2k2+k3T1

The solutions for Ass, Bss, and Css derived from the variable elimination technique still depend on Dss. If we plug in the solutions for the cut species, we can obtain polynomial equations for the remaining species (in this case Dss), but closed form expressions are not necessarily obtainable. In all the cases analyzed in this paper, variables that appear in the analytical solutions for the cut set happen to be approximately constant across a wide range of input values, as they are present in excess relative to other species.

Finally, each parameter group is physically meaningful. For instance, k2k3, k1k3[D]ss2, and k1k2[D]ss2 represent the un-normalized fraction of T1 that exists as A, B, and C, respectively. The normalization factor for these fractions is c/T1, or in this case, simply the sum of all parameter groups. This provides an intuitive way of analyzing how parameter groups affect the overall distribution of T1. For instance, increasing the value of k1 will increase the amount of T1 that exists as B and C, while necessarily decreasing the amount of A (assuming Dss does not change significantly).

2. Wnt model

We analyzed a mathematical model of the canonical Wnt pathway built by Lee et al., 2003. The model is illustrated in Figure 2A, and consists of 7 ODEs and 22 parameters, reproduced in Appendix 1—table 1.

Solving the Wnt model at steady-state

We previously derived an expression for β-catenin in steady-state (Goentoro and Kirschner, 2009):

(A12) [βcat]ss=K171-γ+αuWnt21+4γ1-γ+αuWnt2-1

where the parameters are dimensionless groups of the binding rate constants and protein concentrations:

(A13) α=k4k6k9v14GSK3totAPCtotk5k_6K7K8k13k15
(A14) γ=v12k13K17
(A15) u(Wnt)=1+k3Dvltotk_6k1Wntk2+k1Wnt

The input function u=u(Wnt) corresponds to the rate at which Wnt stimulation inhibits the destruction complex, normalized by k-6. The value of Wnt ranges from 0 to 1 in the model. Please refer to Goentoro and Kirschner (2009) for the physical intuition of each parameter group.

Derivation of linear behavior

We calculate the value of the parameter groups, as well as the value of the input function at saturating Wnt stimulation:

α=66
γ=1.4
uWnt=1=6.0

Within the parameter regime measured in cells, the analytical expression for β-catenin dramatically simplifies. We can perform the following first-order Taylor expansion:

(A16) 1+1+12,1
(A17) ∈=4γ(1γ+αu)2

This holds true for αuγ. Furthermore, we can make the approximation 1-γ+αuαu as long as αu1 also holds. We can encompass these two inequalities within α/u1+γ. The equation simplifies to:

(A18) [βcat]ss K17γ αu

3. ERK model

We analyzed a mathematical model built by Huang and Ferrell (1996), and revised by Sturm et al. (2010). The model is illustrated in Figure 2B, and contains 26 ODEs and 46 parameters, reproduced in Appendix 1—table 2.

We changed two parameters from the original model, which are shown in Appendix 1—table 2. k25 characterizes the negative feedback from dpERK to unphosphorylated Raf, and k27 characterizes the negative feedback from dpERK to phosphorylated Raf. In Sturm et al. (2010), the values of these parameters were estimated, rather than measured. Experimental measurements indicate that dpERK mostly interacts with Raf, and that this feedback causes strong repression of Raf (Dougherty et al., 2005). We therefore increased the value of k25, and set k27 to zero.

Solving the ERK model at steady-state

In the ERK pathway, doubly phosphorylated ERK is produced by the Raf/MEK/ERK cascade of phosphorylation,

(A19) dpERKss=g pRafss

There is a negative feedback within the pathway, such that,

(A20) pRafss=f u,dpERKss

where u is the input function, the concentration of RasGTP (a function of ligand dose).

We first focus on deriving the negative feedback function in Equation A20. Using the variable elimination techniques in section ‘Variable Elimination’, we identify the following cut set:

{Raf, Raf:RasGTP,pRaf,pRaf:P1,MEK:pRaf,pMEK:pRaf,Raf:ppERK,Rafi,Rafi:P4}

This allows us to express the steady-state concentration of pRaf as a function of parameters, and the remaining species in the ERK pathway. Specifically, members of this cut interact directly with, and have dependencies on, the following set:

{P1,MEK,pMEK,dpERK,P4}

With this, we derive the expression for pRafss,

(A21) [pRaf]ss=αuβ[dpERK]ss+(α+γ)u+δRaftot

where the parameter groups are:

(A22) α=k3(k8+kb7)k7[P1]ssk8+
(A23) β=k25(k30+kb29+k29[P4]ss)k29[P4]ssk30+
(A24) γ=k3(k8+kb7)(k9[MEK]ss)k7[P1]ssk8k10+
(A25) δ=k26+kb25k26+

The ellipses indicate additional small terms (i.e. <10% of the previous terms, numerically calculated using the model parameters and u=4.5e4 molecules). All the calculations for this paper use these truncated parameter groups. The complete parameter groups are written below:

α=(k3(k8 + kb7)(k10 + kb9)(k12 + kb11)(k26 + kb25))/([P1]ssk7k8k10k12k26)
β=(k25(k4 + kb3)(k10 + kb9)(k12 + kb11)(k26k30 + k26kb29 +[P4]ssk26k29 +[P4]ssk29k30))/([P4]ssk4k10k12k26k29k30)
γ=(k3(k26 + kb25)(k4k8k10k12 + k4k8k10kb11 + k4k8k12kb9 + k4k10k12kb7 + k4k8kb9kb11 + k4k10kb7kb11 + k4k12kb7kb9 + k4kb7kb9kb11 +[MEK]ssk4k8k9k12 +[MEK]ssk4k8k9kb11 +[MEK]ssk4k9k12kb7 +[MEK]ssk4k9kb7kb11 +[P1]ssk4k7k10k12 +[P1]ssk7k8k10k12 +[P1]ssk4k7k10kb11 +[P1]ssk4k7k12kb9 +[P1]ssk7k8k10kb11 +[P1]ssk7k8k12kb9 +[P1]ssk4k7kb9kb11 +[P1]ssk7k8kb9kb11 + k4k8k10k11[pMEK]ss + k4k8k11kb9[pMEK]ss + k4k10k11kb7[pMEK]ss + k4k11kb7kb9[pMEK]ss))/([P1]ssk4k7k8k10k12k26)
δ= (k4 + kb3k10 + kb9k12 + kb11(k26 + kb25))/(k4k10k12k26)

Physical significance of parameter groups

Next, we would like to develop an intuition for the physical significance of these parameter groups. As discussed in section 1, αu relates to the amount of free, phosphorylated Raf since αu/α+γu+βdpERKss+δ is the fraction of Raf present as pRaf. Thus, as αu increases relative to γu+βdpERKss+δ, the amount of pRaf also increases.

We can define three subpopulations of Raf: Raf inhibited by dpERK, [Ri]; Raf activated by RasGTP (input),[Ra]; and unphosphorylated Raf [Rn]. Specifically:

(A26) Ri=Raf:dpERK+Rafi+Raf:P4
(A27) Ra=pRaf+pRaf:P1+MEK:pRaf+pMEK:pRaf+Raf:RasGTP]
(A28) Rn=[Raf]

We can calculate the steady-state of each subpopulation as:

(A29) [Ri]ss=β[dpERK]ss((α+γ)u+β[dpERK]ss+δ)Raftot
(A30) [Ra]ss=γu((α+γ)u+β[dpERK]ss+δ)Raftot+[pRaf]ss
(A31) [Rn]ss=δ((α+γ)u+β[dpERK]ss+δ)Raftot

Thus, in the same sense that αu relates to the amount of free phosphorylated Raf, βdpERKssrelates to the amount of inhibited Raf, γu relates to the amount of phosphorylated Raf bound to other proteins (not free), and δ relates to the amount of unphosphorylated Raf.

Derivation of linear behavior

Now that we have derived the negative feedback function from Equation A20, we examine Equation A19. The relationship [dpERK]ss=g([pRaf]ss) is analytically intractable, because of the complexity of the phosphorylation cascade. But we know from simulations and experimental observations that it is an ultrasensitive function. From simulations, we find that a 1.3-fold change in pRaf leads to a 9-fold change in dpERK (from 10% to 90% of max, Figure 2—figure supplement 1B-C).

We therefore approximate pRafss by a value Rs within this range, as indicated by the dashed line in Figure 2—figure supplement 1B. Substituting this into the equation above and rearranging, we find that dpERKss becomes a linear function of input:

(A32) [dpERK]ssαβ(RaftotRs1γα)uδβ

Lastly, we write the value of two terms in Equation A32 below, numerically calculated using the parameter values of the model:

αβRaftotRs= 140
αβ1+γα= 13

We can neglect the second term, yielding:

(A33) [dpERK]ssα1βRaftotRsuδβ

Derivation for treating pRaf as a constant

Next, we analyze exactly how the level of pRaf changes with the input u. From earlier, we have that

(A34) dpERKss=gpRafss
(A35) pRaftotss=fu,dpERKss

We can now derive a general expression for the relative change of pRafss with respect to a relative change in u. We use the notation dx^=dlnx=dx/x.

(A36) df^du^=fuuf(1f[dpERK]ssdgd[pRaf]ss)1

Next, we define the response coefficient K between dpERKss andpRafss:

(A37) Kdgd[pRaf]ss[pRaf]ss[dpRaf]ss

From Equation A21, we get the partial derivatives:

(A38) fu=fuβ[dpERK]ss+δβ[dpERK]ss+(α+γ)u+δ
(A39) f[dpERK]ss=fββ[dpERK]ss+(α+γ)u+δ

Using these two equations, we find that:

(A40) df^du^=β[dpERK]ss+δ(1+k)β[dpERK]ss+αu+δ

When K1 and βdpERKss ~ α+γu+δ, we see that

(A41) df^ du^K1

Therefore, pRafss is held constant in the region where the kinase cascade is ultrasensitive and feedback is strong. In this region, it is easy to show that dpERKss becomes a linear function of input.

(A42) dg^g0du^1;       g0= δβ

It is not guaranteed that the system is stable as K increases, but we see from simulations that our parameter regime provides a stable output.

Toy model of the ERK pathway

Here we utilize a toy model to illustrate how ultrasensitivity and strong negative feedback combine to generate input-output linearity. In this model, induction of the output species E is a two-step process:

  1. An input u increases the amount of species R, which in turn influences E as E=gR. There is negative feedback from E to R, which in the limit of strong negative feedback is inversely proportional to E.

  2. Next, we specify the function g(R) such that K=K0, where K is the relative change of E with respect to R. As K0 increases, therefore, the function g(R) becomes more ultrasensitive.

Solving for E, we see that in the limit of K=K01, E becomes a linear function of u, and R is held constant at Rs.

While we do not have an explicit function for g(R) for the full ERK model, we include derivations in section “Derivation for treating pRaf as a constant” that show that these results hold for any function g(R) in the region where K1. We also show that these results hold outside the limit of strong negative feedback, as long as the feedback-inhibited pool of R is comparable to the remaining pool.

Appendix 1—scheme 3
Toy model of the ERK pathway.
https://doi.org/10.7554/eLife.33617.028

4. Tgfβ model

We analyzed a mathematical model built by Schmierer et al. (2008). The model is illustrated in Figure 2C, and consists of 10 ODEs and 14 parameters, reproduced in Appendix 1—tables 3.

Solving the Tgfβ model at steady-state

We use the variable elimination technique described in section ‘Variable Elimination’ to derive an analytical expression for the steady-state concentration of nuclear Smad complex. First, based on the measured parameter values, and as confirmed by simulations, the extent of Smad2-Smad2 binding is limited. We therefore neglect this reaction in subsequent analysis. We identify the following cut of the Tgfβ model:

{S2c,pS2c,S24c,S2n,pS2n,S24n}

which is subject to the conservation equation:

(A43) S2c+pS2c+[S24c]+1a([S2n]+[pS2n]+[S24n])=S2tot

Thus, we can eliminate these variables from the steady-state polynomial solution, with dependence only on variables outside this cut:

{S4c,S4n}

Using this relationship, we derive an expression for the nuclear Smad complex(S24n) at steady-state,

(A44) [S24n]ss=aαu(α+γ)u+βS2tot

where the parameter groups are:

(A45) α=a(kon[S4n]ss+akex2)koff+
(A46) β=PPasekdephoskphosRtotkex2akex2+kin2+
(A47) γ=a(akex2+PPasekdephos)(1akex2+1CIFkin2)+

Here the input function u=u(Tgfβ) is the fraction of receptors activated by Tgfβ ligands. The ellipses indicate additional small terms (i.e. <10% of the previous terms, as calculated using the model parameters, with the variables S4css and S4nss calculated for u=0). All calculations for the paper use these truncated parameter groups. The complete parameter groups are written below:

α=(a([S4n]sskoff + CIF[S4n]sskin2 + CIFPPase[S4c]sskdephos + CIF[S4c]ss[S4n]sskon + CIF[S4c]ssakex2))/(CIF[S4c]sskoff)
β=(PPasekdephos(kin2 + akex2)(koff + CIFkin2 + CIF[S4c]sskon))/(CIFRtot[S4c]sskex2konkphos)
γ=((PPasekdephos + akex2)(kin2koff + CIFkin22 + akex2koff + CIF[S4c]sskin2kon + CIFakex2kin2 +[S4c]ssakex2kon))/(CIF[S4c]sskex2kin2kon)

Physical significance of parameter groups

Next, we would like to develop an intuition for the physical significance of these parameter groups. As discussed in section 1, αu relates to the amount of nuclear Smad complex since αu/α+γu+β is the fraction of Smad2 present as S24n. Thus, as αu increases relative toγu+β, the amount of S24n also increases.

By definition, the parameter groups β and γu capture the remaining input-independent and input-dependent polynomials, respectively. Nevertheless, we would like to understand the physical significance of the parameter groups. We can calculate the amount of unphosphorylated Smad2 as:

(A48) [S2c]ss+1a[S2n]ss=β+δuβ+(α+γ)uS2tot
(A49) δ=PPasekdephoskoff+CIFkin2+CIF[S4c]sskonCIF[S4c]sskex2kon

δ captures the dependence of nuclear, unphosphorylated Smad on the input. With the measured parameters, βδu, so we have

(A50) [S2c]ss+1a[S2n]ssββ+(α+γ)uS2tot

This means that β relates to the amount of unphosphorylated Smad2 in the same sense that αu relates to nuclear Smad complex. We can also express the remaining Smad2 species as:

(A51) [pS2c]ss+[S24c]ss+1a[pS2n]ss=(γδ)u((α+γ)u+β)S2tot

However, as δ is of the same order of magnitude as γ, the parameter group γ only loosely relates to these remaining species of Smad2.

Derivation of linear behavior

Within the parameter values measured in cells, the behavior of Smad complex dramatically simplifies. Using the measured values (Appendix 1—table 3), the parameter groups are

αu=3.1
γu=1.3
β=46

where we have used a non-saturating input (u= 0.2). Therefore, with measured parameters, βα+γu. With this, the denominator in the [S4n]ss equation simplifies, and the concentration of Smad complex becomes a linear function of the input:

(A52) [S24n]ssαS2totβu
Appendix 1—table 1
Parameters, variables, and equations of the Wnt model.
https://doi.org/10.7554/eLife.33617.029
ParameterLabelValue
Activation rate of Disheveled/Dvl by Wnt

k1

0.182

min-1

Inactivation rate of Dvl

k2

1.82102

min-1

Dissociation of destruction complex (DC) by active Dvl

k3

5.00102

nM-1 min-1

Phosphorylation of DC

k4

0.267

min-1

Dephosphorylation of DC

k5

0.133

min-1

Forward rate for DC binding

k6

9.09102

nM-1 min-1

Reverse rate for DC binding

k-6

0.909

min-1

Dissociation constant for APC:axin binding

K7

50

nM

Dissociation constant for β-catenin:DC binding

K8

120

nM

Phosphorylation rate of β-catenin

k9

206

min-1

Rate of phosphorylated β-catenin release from DC

k10

206

min-1

Degradation rate of phosphorylated β-catenin

k11

0.417

min-1

Synthesis rate of β-catenin

v12

0.423

nM min-1

Degradation rate of β-catenin

k13

2.57104

min-1

Synthesis rate of axin

v14

8.22105

nM min-1

Degradation rate of axin

k15

0.167

min-1

Dissociation constant for β-catenin:TCF binding

K16

30

nM

Dissociation constant for β-catenin:APC binding

K17

1200

nM

Total concentration of Disheveled

Dvltot

100

nM

Total concentration of adenomatous polyposis coli

APCtot

100

nM

Total concentration of T-cell factor

TCFtot

15

nM

Total concentration of glycogen synthase kinase 3β

GSK3tot

50

nM

Independent VariableLabel
Active Disheveled

X2

APC*/axin*/GSK3 (* denotes phosphorylated)

X3

APC/axin/GSK3

X4

β-catenin*/APC*/axin*/GSK3

X9

β-catenin*

X10

β-catenin

X11 (βcat)

axin

X12

Dependent VariableLabel
Inactive Disheveled

X1

GSK3

X5

APC/axin

X6

APC

X7

β-catenin/APC*/axin*/GSK3

X8

TCF

X13

β-catenin/TCF

X14

β-catenin/APC

X15

Differential Equations

[X2]˙=k1Wnt(Dvltot[X2])k2[X2]

(1+[X11]K8)[X3]˙+[X3]K8[X11]˙=k4[X4]k5[X3]k9[X3][X11]K8+k10[X9]

[X4]˙=(k3[X2]+k4+k6)[X4]+k5[X3]+k6GSK3totK17[X12]APCtotK7(K17+[X11])

[X9]˙=k9[X3][X11]K8k10[X9]

[X10]˙=k10[X9]k11[X10]

(1+[X3]K8+K16TCFtot(K16+[X11])2+K17APCtot(K17+[X11])2)[X11]˙+[X11]K8[X3]˙=v12(k9[X3]K8+k13)[X11]



(1+K17APCtotK7(K17+[X11]))[X12]˙K17[X12]APCtotK7(K17+[X11])2[X11]˙=k3[X2][X4]k6GSKtotK17[X12APCtot]K7(K17+[X11])+k6[X4]+v14k15[X12]


Equations for fast equilibrium reactions

X1=Dvltot-X2

X5=GSK3tot


[X6]=K17[X12APCtot]K7(K17+[X11])

[X7]=K17APCtotK17(K17+[X11])

[X8]=[X3][X11]K8

[X13]=K16TCFtotK16+[X11]

[X14]=[X11]TCFtotK16+[X11]

X15=X11APCtotK17+X11
Appendix 1—table 2
Parameters, variables, and equations of the ERK model.

Values highlighted in yellow have been changed from the original model (explained in section ‘ERK Model’).

https://doi.org/10.7554/eLife.33617.030
ParameterLabelValue
Forward rate for Raf:RasGTP binding

k3

1.67106

molecule1s1

Reverse rate for Raf:RasGTP binding

kb3

5.3103

s-1

Phosphorylation rate for Raf by RasGTP

k4

1

s-1

Forward rate of pRaf:P1 binding

k7

1.18104

molecule-1s-1

Reverse rate of pRaf:P1 binding

kb7

0.2

s-1

Dephosphorylation rate of pRaf by P1

k8

1

s-1

Forward rate of MEK:pRaf binding

k9

1.95105

molecule-1s-1

Reverse rate of MEK:pRaf binding

kb9

3.3102

s-1

Phosphorylation rate of MEK by pRaf

k10

3.5

s-1

Forward rate of pMEK:pRaf binding

k11

1.95105

molecule-1s-1

Reverse rate of pMEK:pRaf binding

kb11

3.3102

s-1

Phosphorylation rate of pMEK by pRaf

k12

2.9

s-1

Forward rate of dpMEK:P2 binding

k13

2.38105

molecule-1s-1

Reverse rate of dpMEK:P2 binding

kb13

0.8

s-1

Dephosphorylation rate of dpMEK by P2

k14

5.8102

s-1

Forward rate of pMEK:P2 binding

k15

4.5107

molecule-1s-1

Reverse rate of pMEK:P2 binding

kb15

0.5

s-1

Dephosphorylation rate of pMEK by P2

k16

5.8102

s-1

Forward rate of ERK:dpMEK binding

k17

8.9105

molecule-1s-1

Reverse rate of ERK:dpMEK binding

kb17

1.83102

s-1

Phosphorylation rate of ERK by dpMEK

k18

16

s-1

Forward rate of pERK:dpMEK binding

k19

8.9105

molecule-1s-1

Reverse rate of pERK:dpMEK binding

kb19

1.83102

s-1

Phosphorylation rate of pERK by dpMEK

k20

5.7

s-1

Forward rate of pERK:P3 binding

k21

8.33106

molecule-1s-1

Reverse rate of pERK:P3 binding

kb21

0.5

s-1

Dephosphorylation rate of pERK by P3

k22

0.246

s-1

Forward rate of dpERK:P3 binding

k23

2.35105

molecule-1s-1

Reverse rate of dpERK:P3 binding

kb23

0.6

s-1

Dephosphorylation rate of dpERK by P3

k24

0.246

s-1

Forward rate of Raf:dpERK binding

k25

1106

molecule-1s-1

Reverse rate of Raf:dpERK binding

kb25

1

s-1

Hyper-phosphorylation rate of Raf by ppERK

k26

10

s-1

Forward rate of pRaf:dpERK binding

k27

0

molecule-1s-1

Reverse rate of pRaf:dpERK binding

kb27

1

s-1

Hyper-phosphorylation rate of phosphorylated Raf by dpERK

k28

10

s-1

Forward rate of Rafi:P4 binding

k29

5105

molecule-1s-1

Reverse rate of Rafi:P4 binding

kb29

0.2

s-1

Dephosphorylation rate of Rafi by P4

k30

0.5

s-1

Total Raf

Raftot

4104

molecules

Total MEK

MEKtot

2.1107

molecules

Total ERK

ERKtot

2.21107

molecules

Total phosphatase P1

P1tot

4104

molecules

Total phosphatase P2

P2tot

4105

molecules

Total phosphatase P3

P3tot

1107

molecules

Total phosphatase P4

P4tot

4104

molecules

VariableLabel
Unphosphorylated Raf

Raf

Raf bound to RasGTP

Raf:RasGTP

Phosphorylated Raf

pRaf

Phosphatase for phosphorylated Raf

P1

Phosphorylated Raf bound to its phosphatase

pRaf:P1

Unphosphorylated MEK

MEK

MEK bound to its kinase

MEK:pRaf

Phosphorylated MEK

pMEK

Phosphorylated MEK bound to its kinase

pMEK:pRaf

Doubly-phosphorylated MEK

dpMEK

MEK phosphatase

P2

Doubly-phosphorylated MEK bound to its phosphatase

dpMEK:P2

Phosphorylated MEK bound to its phosphatase

pMEK:P2

Unphosphorylated ERK

ERK

ERK bound to its kinase

ERK:dpMEK

Phosphorylated ERK

pERK

Phosphorylated ERK bound to its kinase

pERK:dpMEK

Doubly-phosphorylated ERK

dpERK

ERK phosphatase

P3

Phosphorylated ERK bound to its phosphatase

pERK:P3

Doubly-phosphorylated ERK bound to its phosphatase

dpERK:P3

Raf bound to doubly-phosphorylated ERK

Raf:dpERK

Hyper-phosphorylated, ‘inactive’ Raf

Rafi

Phosphorylated Raf bound to doubly-phosphorylated ERK

pRaf:dpERK

Phosphatase for hyper-phosphorylated Raf

P4

Hyper-phosphorylated Raf bound to its phosphatase

Rafi:P4

Differential Equations

[Raf]˙=k3[Raf]u(EGF)+kb3[Raf:RasGTP]+k8[pRaf:P1]k25[Raf][dpERK]+kb25[Raf:dpERK]+k30[Rafi:P4]

[Raf:RasGPT]˙=k3[Raf]u(EGF)(kb3+k4)[Raf:RasGTP]

[pRaf]˙=k4[Raf:RasGTP]k7[pRaf][P1]+kb7[pRaf:P1]k9[MEK][pRaf]+(kb9+k10)[MEK:pRaf]k11[pMEK][pRaf]+(kb11+k12)[pMEK:pRaf]k27[pRaf][dpERK]+kb27[pRaf:dpERK]

[P1]˙=k7[pRaf][P1]+(kb7+k8)[pRaf:P1]

[pRaf:P1]˙=k7[pRaf][P1](kb7+k8)[pRaf:P1]

[MEK]˙=k9[MEK][pRaf]+kb9[MEK:pRaf]+k16[pMEK:P2]

[MEK:pRaf]˙=k9[MEK][pRaf](kb9+k10)[MEK:pRaf]

[pMEK]˙=k10[MEK:pRaf]k11[pMEK][pRaf]+kb11[pMEK:pRaf]+k14[dpMEK:P2]k15[pMEK][P2]+kb15[pMEK:P2]

[pMEK:pRaf]˙=k11[pMEK][pRaf](kb11+k12)[pMEK:pRaf]

[dpMEK]˙=k12[pMEK:pRaf]k13[dpMEK][P2]+kb13[dpMEK:P2]k17[ERK][dpMEK]+(kb17+k18)[ERK:dpMEK]k19[pERK][dpMEK]+(kb19+k20)[pERK:dpMEK]

[P2]˙=k13[dpMEK][P2]+(kb13+k14)[dpMEK:P2]k15[pMEK][P2]+(kb15+k16)[pMEK:P2]

[dpMEK:P2]˙=k13[dpMEK][P2](kb13+k14)[dpMEK:P2]

[pMEK:P2]˙=k15[pMEK][P2](kb15+k16)[pMEK:P2]

[ERK]˙=k17[ERK][dpMEK]+kb17[ERK:dpMEK]+k22[pERK:P3]

[ERK:dpMEK]˙=k17[ERK][dpMEK](kb17+k18)[ERK:dpMEK]

[pERK]˙=k18[ERK:dpMEK]k19[pERK][dpMEK]+kb19[pERK:dpMEK]k21[pERK][P3]+kb21[pERK:P3]+k24[dpERK:P3]

[pERK:dpMEK]˙=k19[pERK][dpMEK](kb19+k20)[pERK:dpMEK]

[dpERK]˙=k20[pERK:dpMEK]k23[dpERK][P3]+kb23[dpERK:P3]k25[Raf][dpERK]+(kb25+k26)[Raf:dpERK]k27[pRaf][dpERK]+(kb27+k28)[pRaf:dpERK]

[P3]˙=k21[pERK][P3]+(kb21+k22)[pERK:P3]k23[dpERK][P3]+(kb23+k24)[dpERK:P3]

[pERK:P3]˙=k21[pERK][P3](kb21+k22)[pERK:P3]

[dpERK:P3]˙=k23[dpERK][P3](kb23+k24)[dpERK:P3]

[Raf:dpERK]˙=k25[Raf][dpERK](kb25+k26)[Raf:dpERK]

[Rafi]˙=k26[Raf:dpERK]+k28[pRaf:dpERK]k29[Rafi][P4]+kb29[Rafi:P4]

[pRaf:dpERK]˙=k27[pRaf][dpERK](kb27+k28)[pRaf:dpERK]

[P4]˙=k29[Rafi][P4]+(kb29+k30)[Rafi:P4]

[Rafi:P4]˙=k29[Rafi][P4](kb29+k30)[Rafi:P4]
Algebraic Equations for conserved species

Raftot=[Raf]+[Raf:RasGTP]+[pRaf]+[pRaf:P1]+[MEK:pRaf]+[pMEK:pRaf]+[Raf:dpERK]+[Rafi]+[pRaf:dpERK]+[Rafi:P4]

MEKtot=[MEK]+[MEK:pRaf]+[pMEK]+[pMEK:pRaf]+[dpMEK]+[dpMEK:P2]+[pMEK:P2]+[ERK:dpMEK]+[pERK:dpMEK]

ERKtot=[ERK]+[ERK:dpMEK]+[pERK]+[pERK:dpMEK]+[dpERK]+[pERK:P3]+[dpERK:P3]+[Raf:dpERK]+[pRaf:dpERK]

P1tot=[P1]+[pRaf:P1]

P2tot=[P2]+[dpMEK:P2]+[pMEK:P2]

P3tot=[P3]+[pERK:P3]+[dpERK:P3]

P4tot=[P4]+[Rafi:P4]
Appendix 1—table 3
Parameters, variables and equations of the Tgfβ model.
https://doi.org/10.7554/eLife.33617.031
ParameterLabelValue
Phosphorylation rate of Smad2

kphos

4.0104

nM1s1

Dephosphorylation rate of Smad2

kdephos

6.6103

nM1s1

Nuclear import rate of Smad2

kin2

2.6103

s1

Nuclear export rate of Smad2

kex2

5.6103

s1

Nuclear import rate of Smad4

kin4

2.6103

s1

Nuclear export rate of Smad4

kex4

2.6103

s1

Smad complex import factor

CIF

5.7

Forward rate for Smad complex binding

kon

1.8103

nM1s1

Reverse rate for Smad complex binding

koff

1.6102

s1

Cytoplasmic to nuclear volume ratio

a

2.3

Total Smad2 (initialized to cytoplasm)

S2tot

73.0

nM

Total Smad4 (initialized to cytoplasm)

S4tot

73.0

nM

Total phosphatase in nucleus

PPase

1

nM

Total Receptors

Rtot

1

nM

VariableLabel
Cytoplasmic Smad2

S2c

Cytoplasmic phosphorylated Smad2

pS2c

Cytoplasmic Smad4

S4c

Cytoplasmic Smad2:Smad4 complex

S24c

Cytoplasmic Smad2:Smad2 complex

S22c

Nuclear Smad2

S2n

Nuclear phosphorylated Smad2

pS2n

Nuclear Smad4

S4n

Nuclear Smad2:Smad4 complex

S24n

Nuclear Smad2:Smad2 complex

S22n

Differential Equations

[S2c]˙=kphosu(Tgfβ)[S2c]kin2[S2c]+kex2[S2n]

[pS2c]˙=kphosu(Tfgβ)[S2c]kin2[pS2c]kon[pS2c]([S4c]+2[pS2c])+koff([S24c]+2[S22c])+kex2[pS2n]

[S4c]˙=kin4[S4c]kon[pS2c][S4c]+koff[S24c]+kex4[S4n]

[S24c]˙=kon[pS2c][S4c]koff[S24c]kin2CIF[S24c]

[S22c]˙=kon[pS2c]2koff[S22c]kin2CIF[S22c]

[S2n]˙=akin2[S2c]akex2[S2n]+kdephosPPase[pS2n]

[pS2n]˙=akin2[pS2c]akex2[pS2n]kdephosPPase[pS2n]kon[pS2n]([S4n]+2[pS2n])+koff([S24n]+2[S22n])

[S4n]˙=akin4[S4c]akex4[S4n]kon[pS2n][S4n]+koff[S24n]

[S24n]˙=akin2CIF[S24c]+kon[pS2n][S4n]koff[S24n]

[S22n]˙=akin2CIF[S22c]+kon[pS2n]2koff[S22n]
Algebraic Equations for conserved species

S2tot=[S2c]+[pS2c]+[S24c]+2[S22c]+(2[S22n]+[S24n]+[pS2n]+[S2n])

S4tot=[S4c]+[S24c]+1a([S24n]+[S4n])
Appendix 1—table 4
Examples of biological systems where the Wnt, ERK, and Tgfβ pathways have been shown to produce graded response in single-cell level.
https://doi.org/10.7554/eLife.33617.032
PathwaySystems where graded response has been observedReferences
Live imaging of single cells
Tgfβ pathwayMouse myoblastsFrick et al., 2017; Warmflash et al. (2012)
Human epidermal keratinocytesNicolás et al. (2004); Warmflash et al. (2012) ; Schmierer et al. (2008); Vizán et al. (2013)
Human cervical epithelial cellsNicolás et al. (2004)
Human breast epithelial cellsStrasen et al. (2018)
Canonical Wnt pathwayHuman embryonic kidney cellsKafri et al. (2016) (this is the only published live single-cell imaging study in the Wnt pathway so far)
ERK pathwayMouse fibroblastsToettcher et al. (2013)
Mouse embryonic fibroblastsMackeigan et al. (2005)
Human non-small cell lung carcinomaCheong et al., 2011
Human mammary gland cellsSelimkhanov et al. (2014); Perrett et al., 2013Perrett et al., 2013
Human cervical epithelial cellsVoliotis et al. (2014); Whitehurst et al. (2004); Perrett et al., 2013Perrett et al., 2013
Human foreskin fibroblastsWhitehurst et al. (2004)
Immunofluorescence and FACS studies
Tgfβ pathwayXenopus embryoSchohl and Fagotto (2002)
Mouse testesItman et al. (2009)
Zebrafish embryoDubrulle et al. (2015)
Canonical Wnt pathwayXenopus embryoSchneider et al. (1996); Fagotto and Gumbiner (1994); Schohl and Fagotto (2002)
Mouse embyoAulehla et al., 2008
PlanariaSureda-Gómez et al., 2016
Sea anemone embryoWikramanayake et al., 2003
ERK pathwayChick embryoDelfini et al. (2005)
Xenopus embryoSchohl and Fagotto (2002)
Human T lymphocyte cellsLin et al. (2009)
Rat adrenal gland cellsSantos et al., 2007

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
    Beta-catenin localization during Xenopus embryogenesis: accumulation at tissue and somite boundaries
    1. F Fagotto
    2. BM Gumbiner
    (1994)
    Development 120:3667–3679.
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
    Relationship between epidermal growth factor receptor occupancy and mitogenic response. quantitative analysis using a steady state model system.
    1. DJ Knauer
    2. HS Wiley
    3. DD Cunningham
    (1984)
    The Journal of Biological Chemistry 259:5623–5654.
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
    Allosteric proteins as logarithmic sensors
    1. N Olsman
    2. L Goentoro
    (2016)
    Proceedings of the National Academy of Sciences 113:E4423–E4430.
    https://doi.org/10.1073/pnas.1601791113
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
    β-catenin, MAPK and Smad signaling during early Xenopus development
    1. A Schohl
    2. F Fagotto
    (2002)
    Development 129:37–52.
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70
  71. 71
  72. 72
  73. 73
  74. 74
  75. 75
  76. 76
  77. 77
  78. 78
  79. 79
  80. 80
  81. 81

Decision letter

  1. Wenying Shou
    Reviewing Editor; Fred Hutchinson Cancer Research Center, United States
  2. Aviv Regev
    Senior Editor; Broad Institute of MIT and Harvard, United States
  3. Wenying Shou
    Reviewer; Fred Hutchinson Cancer Research Center, United States
  4. Steven S Andrews
    Reviewer; Fred Hutchinson Cancer Research Center, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your work entitled "Signaling pathways as linear transmitters" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Wenying Shou as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by a Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Steven S Andrews (Reviewer #2).

Our decision has been reached after consultation between the reviewers. We all found your work interesting, but we all had problems with the precise meaning of "physiological ranges".

Reviewer #1:

Nunns and Goentoro examined the conserved cores of three signaling pathways (Wnt, ERK, and TGFb). Using mathematical modeling, they found that despite differences in pathway architecture, the three pathways all behaved like linear signal transmitters within physiological ranges. Using experiments to test two pathways (Wnt and ERK), they indeed found linear input-output relationship.

I like how they used dimensionless analysis to gain insights on mechanisms of linearity. My major comments are:

1) They need to describe more carefully the "physiology range": what do they mean by physiological? Moreover, I assume that every parameter has its own range, and dimensionless numbers comprising many parameters probably have a quite large range. What are their ranges? Does most of the range fall within "physiological range"?

2) On a related note, for ERK which is known to sometimes exhibit switch-like responses, how would those parameters and inputs fare in your analytics? Would they have told you that the response will be non-linear?

3) Is there any way to make the wiring diagrams more intuitive and the important points more obvious? Perhaps color coding those parameters in α or γ? More details in legend? For example, why α represents β-catenin degradation by DC is not clear: k11 – which I assume means degradation of phosphor β catenin – is not in α. Why not? The information might already be in supplementary data, but intuition about an expression should be in the main text if you want to appeal to a broad audience.

Reviewer #2:

The authors used modeling and then experiment to show linear input-output relations in several model signaling pathways, which were the canonical Wnt, ERK, and Tgf-β pathways. The modeling started with established ODE-based models and reduced them to simplified analytical expressions that related output to input; these expressions did not show linear relationships in general but did in physiologically relevant parameter regimes. The experiments showed linearity in the Wnt and ERK pathways, and that this linearity was reduced through system perturbations.

This is very good work on an important problem. The methods and analysis were appropriate, and the paper is well written.

1) It should be clarified whether the authors are referring to linear signal transmission relative to the signaling pathway ligand concentrations, or to the fraction of bound receptors. This issue arises because the authors show linearity relative to the ligand concentrations in the model analysis portion of this work, but then linearity relative to the fraction of bound receptors in the experimental work on the Wnt pathway. More specifically, it appears that the use of Wnt as the x-axis in Figure 2D contradicts statements in subsection “Linearity in the Wnt and ERK pathways was observed experimentally”, which say that "linearity does not extend upstream to Wnt dose".

2) My experience with cell signaling pathway ODE models is that they capture the known biological interactions reasonably correctly, they agree well with a specific set of test data, and they get the signaling dynamics qualitatively correct, but they are rarely accurate enough for quantitative predictions on untested problems. However, this work uses three models for quantitative predictions. On the one hand, I applaud this approach because it fulfills a major objective of modeling research. On the other hand, I feel that more work is required to convince the reader that the models are in fact accurate enough to warrant the conclusions that are reached from them. This includes the models having sufficient accuracy over both the parameter and time ranges that are considered. For example, if the grey lines in Figure 3 were model calculations (without any additional fitting), I would have greater trust in the models.

3) A persistent concern that I had is that most dose-response functions for cell signaling can be modeled well by Hill functions, and Hill functions are linear relationships in the small dose regime. Is the linearity observed here simply an observation of Hill functions in the small dose regime? If so, then all of these results are reasonably trivial. As part of addressing this issue, it would help if the authors could give the physiologically typical concentration ranges for Wnt, EGF, and Tgf-β.

Reviewer #3:

This manuscript aims to test whether signaling pathways have linear input-output response characteristics. A simplification of existing mathematical models indicates this being the case for the canonical Wnt, ERK and TGF-β pathways, despite striking differences in core network architectures among these pathways. It is interesting to see how different complex pathways can be approximated by linear systems. The manuscript concludes by experiments (Western Blots) testing the theoretical results, as well as the breakdown of linearity as some key features of the pathways are altered.

In summary the major point here seems to be that in many of these diverse, complicated, pathways, a linear dose-response can be observed as an approximation, and this might be borne out in experimental conditions as well. In fact, not only can the linear dose response be observed, it appears rather robust to perturbations until you exceed some critical value. While these results are potentially interesting, there are major weaknesses that should be addressed before the manuscript could be considered for publication in eLife.

1) Besides the three pathways chosen, there are other pathways that are similarly well studied, both experimentally as well as mathematically. The NF-kappaB pathway is a good example. Why was the NF-kappaB pathway or any other pathway not considered? You could make a statement that proving the generality of this finding will require investigating additional pathways in the future.

2) A linear input-output relationship due to negative feedback has already been described for a MAP Kinase pathway, see PMID:19079053. How are the current results novel compared to these earlier findings (besides studying a human MAPK pathway)?

3) The claims of linearity in this manuscript are mainly based on visual examination. However, there are rigorous ways to measure linearity, based on the L1-norm, as described in PMID:19279212 and PMID:23385595, which should be cited. It would be necessary to apply the L1-norm throughout all figures of the manuscript to make computational and experimental claims of linearity quantitative.

4) For deriving formula [4], apha/u>>1 is necessary. However, based on the SI, α/u=11, and it is actually smaller for some of the range plotted in Figure 2. Then the linearity arises from 'α / u >> 1 + γ'. One could/should explore when α = (1+γ)/u * S where S is a scaling factor set to say, 0.9, 1, 1.1, 2, 10, and 100. At 10 and 100, one might expect the condition to hold, but at S around 1 the simplification probably fails. Do we then see lack of linearity? One can then re-arrange the terms for each of the other terms (u or γ) and explore similarly. The same argument can be made for the Erk and Tgf[β] pathways.

4) How robust is linearity? To answer this question, two modifications are needed: (i) extend the range of u for each plot in Figure 2, showing the nonlinearity of the curves, then indicate the linear range on such nonlinear curves, e.g. with a different color; (ii) Change a parameter and plot a family of curves, indicating the linear range (determined based on the L1-norm) on each individual curve. This can be done for a couple of parameters, but it is especially important for the parameters mediating interactions that cause linearity to break down (described in subsection “Linearity in the Wnt and ERK pathways is modulated by perturbation to parameters”).

5) If linear responses indeed help cellular signal processing as a result of convergent evolution then linearity should probably occur in single cells that process inputs. Unfortunately, the Western Blots in Figure 3 only test linearity at the population average level. The linearity of the average does not imply linearity in single cells. So, are single cell input-output response characteristics linear? This should be checked at least for one signaling pathway, possibly by new experiments or otherwise using data by others, see for example PMID:25504722. In fact, recent papers indicate that some pathways' responses are dynamic/oscillatory, noisy or bimodal at the single cell level. All of this should be addressed/discussed. Where is then the linearity?

6) What defines a "physiological range"? As noted, ultrasensitivity has been described for the ERK pathway's input response. Ultrasensitivity seems contradictory to linearity. Is then ultrasensitivity not physiological? Moreover, what is physiological in a given tissue or organism (from yeast to frog eggs to worms to humans) might not be physiological in another one. All of this should be discussed and clarified.

7) Considering the importance of fold-change detection in biology, it would be important to assess the number of decades over which linearity holds for each pathway. In fact, any nonlinear function can be approximated by a linear relationship except near extremum points (this is the essence of the Taylor approximation). Are these observations more than Taylor expansion? Maybe we can tell based on the decades over which linearity holds.

8) The data in Figure 3 should be expanded into the domains of nonlinearity/saturation and then the range of linearity should be indicated on top of such curves (linearity measured using the L1-norm). This should be done for all panels. In fact, without showing the full curve (including the nonlinear parts), panel 3C does not indicate that linearity is lost or that its range shrinks.

9) What ranges on the theoretical plots do the experimental results correspond to? Theory and experiment should be better connected.

10) Experimental verification is not trivial in these systems. First, the "quantitative' Western blots are not truly quantitative, despite ensuring that fluorescent antibodies are within linear range, etc. Western blots, by their very nature, are notoriously noisy, and therefore require a lot of "normalization" which might artifactually introduce the appearance of linearity. If, for example, the denominator is very large for normalization, the signal can appear more linear. Say, if you compare y=x and y=x^2, but then divide by a large factor K such that within the range of 'x' being explored K>>y, then both equations will be nearly 0 and will likely show a "linear" response (also true according to the Taylor expansion near 0).

11) The variability in the Western blots is borne out in Figure 3A where the various 'Wnt' doses have very variable pLRP5/6 response (input). Likewise, in Figure 3B, for EGF. These few points do not necessarily demonstrate linearity other than the fact that the linear approximation can be applied to just about any relationship. In addition, comparing Figure 3A with Figure 3D where linearity breaks down with a Raf mutation preventing feedback, it still appears linear until 3ng/mL, a dose that is not reported for the wild-type (Figure 3A). In fact, examining both plots between 0 to 2ng/mL, one could conclude instead that both systems are 'linear'.

12) Same with the Wnt pathway. Here it is already recognized that linearity was not lost, but this is attributed to side effects of the drug used – something that could be tested perhaps with a different drug or a mutation similar to the Erk pathway analysis.

https://doi.org/10.7554/eLife.33617.034

Author response

We have now addressed the editors’ and reviewers’ comments, and believe that the manuscript has been made stronger and more precise by the process. We therefore would love to have the manuscript reconsidered, as we believe eLife would provide a good home for reaching the suitable audience for the work. Specifically, we have addressed the primary concern regarding the precise meaning of “physiological range”. We agree with the editor and reviewers that this gives a misleading impression that we have information about the parameter ranges across all biological contexts. We have now revised the misleading phrase into “measured parameters in (the name of the specific systems)”, “physiologically relevant parameter ranges”, or “some physiological contexts”.

Reviewer #1:

Nunns and Goentoro examined the conserved cores of three signaling pathways (Wnt, ERK, and TGFb). Using mathematical modeling, they found that despite differences in pathway architecture, the three pathways all behaved like linear signal transmitters within physiological ranges. Using experiments to test two pathways (Wnt and ERK), they indeed found linear input-output relationship.

I like how they used dimensionless analysis to gain insights on mechanisms of linearity. My major comments are:

1) They need to describe more carefully the "physiology range": what do they mean by physiological? Moreover, I assume that every parameter has its own range, and dimensionless numbers comprising many parameters probably have a quite large range. What are their ranges? Does most of the range fall within "physiological range"?

We see now that the phrasing “physiological range” gives a misleading impression that we have information about the physiological range of the parameters across all biological contexts. To avoid this confusion, we now described it more clearly in the manuscript that we considered how the analytical solutions behaved in a physiologically relevant parameter regime – i.e., the values of parameters that have been measured or estimated in specific biological contexts. All changes are in marked in red. As linearity emerges with these measured parameters, it argues for the physiological relevance of the predicted linearity. Additionally, we experimentally observed linearity in the two human cell lines.

While we have the parameter values that have been measured in some biological contexts, we do not have the full range of the parameter groups across physiological contexts. However, as also requested by reviewer 3 (Point 5), we now include robustness analysis that shows that linearity occurs through a considerable parameter range in the models, although not limitless. It therefore further argues the significance that the measured parameters fall within the range where linearity prevails.

2) On a related note, for ERK which is known to sometimes exhibit switch-like responses, how would those parameters and inputs fare in your analytics? Would they have told you that the response will be non-linear?

Switch-like response in the ERK model could be achieved by decreasing β or increasing 𝛼𝛼, both amount to weakening the negative feedback. The analysis would tell us that the response is nonlinear (see Figure 2—figure supplement 6C), i.e., L-1 norm > 0.1. Consistently, we observed experimentally that mutating Raf that mediates the negative feedback converts linear to nonlinear response (Figure 3D).

3) Is there any way to make the wiring diagrams more intuitive and the important points more obvious? Perhaps color coding those parameters in α or γ? More details in legend? For example, why α represents β-catenin degradation by DC is not clear: k11 – which I assume means degradation of phosphor β catenin – is not in α. Why not? The information might already be in supplementary data, but intuition about an expression should be in the main text if you want to appeal to a broad audience.

We thank the reviewer for this suggestion. We have revised Figure 2. We added in the main text intuitive description of the parameter grouping, with references to the shading. We also revised the legend to explain why certain parameters are not in the parameter groups (e.g., k11 describes an irreversible degradation of phosphorylated β-catenin, and therefore does not affect the level of unphosphorylated β-catenin we are solving for).

Reviewer #2:

The authors used modeling and then experiment to show linear input-output relations in several model signaling pathways, which were the canonical Wnt, ERK, and Tgf-β pathways. The modeling started with established ODE-based models and reduced them to simplified analytical expressions that related output to input; these expressions did not show linear relationships in general but did in physiologically relevant parameter regimes. The experiments showed linearity in the Wnt and ERK pathways, and that this linearity was reduced through system perturbations.

This is very good work on an important problem. The methods and analysis were appropriate, and the paper is well written.

1) It should be clarified whether the authors are referring to linear signal transmission relative to the signaling pathway ligand concentrations, or to the fraction of bound receptors. This issue arises because the authors show linearity relative to the ligand concentrations in the model analysis portion of this work, but then linearity relative to the fraction of bound receptors in the experimental work on the Wnt pathway. More specifically, it appears that the use of Wnt as the x-axis in Figure 2D contradicts statements in subsection “Linearity in the Wnt and ERK pathways was observed experimentally”, which say that "linearity does not extend upstream to Wnt dose".

We agree this could be confusing since the form of the input function u(ligand) is different across the models, and not necessarily equal to ligand concentration. In the Wnt pathway, u(ligand) describes the Wnt-dependent action of Dvl in promoting dissociation of the destruction complex. In the Tgfβ pathway, u(ligand) is the fraction of active receptors. In the ERK pathway, u(ligand) is the concentration of EGF-activated Ras-GTP. To clarify this, we have now made u(ligand) more obvious in two ways:

A) We include u(ligand) function in the x-axis of Figure 2D-F and add brief intuitive description of each u(ligand), as opposed to having this only in the legend previously.

B) We include in the wiring diagrams in Figure 2A-C illustration of the upstream processes that constitute the function u(ligand), to further emphasize that they are not simply ligand concentration.

2) My experience with cell signaling pathway ODE models is that they capture the known biological interactions reasonably correctly, they agree well with a specific set of test data, and they get the signaling dynamics qualitatively correct, but they are rarely accurate enough for quantitative predictions on untested problems. However, this work uses three models for quantitative predictions. On the one hand, I applaud this approach because it fulfills a major objective of modeling research. On the other hand, I feel that more work is required to convince the reader that the models are in fact accurate enough to warrant the conclusions that are reached from them. This includes the models having sufficient accuracy over both the parameter and time ranges that are considered. For example, if the grey lines in Figure 3 were model calculations (without any additional fitting), I would have greater trust in the models.

We agree that we need to clarify that we are not simply observing linearity in a small dose regime. We now explain more explicitly in the text and include Figure 2—figure supplement 2 to emphasize that the predicted linearity extends the entire or almost the entire dynamic rangeof the systems.

3) A persistent concern that I had is that most dose-response functions for cell signaling can be modeled well by Hill functions, and Hill functions are linear relationships in the small dose regime. Is the linearity observed here simply an observation of Hill functions in the small dose regime? If so, then all of these results are reasonably trivial. As part of addressing this issue, it would help if the authors could give the physiologically typical concentration ranges for Wnt, EGF, and Tgf-β.

Based on our own experience, we hold the same view that when working across biological contexts, we use the models of signaling pathways for exploring qualitative, systemslevel behaviors. With regard to this, the three models have track record of success:

- The ERK model predicts ultrasensitive behavior (Huang and Ferrell, 1996).

- The Wnt model enables deducing the differential roles of the two scaffolds in the pathway (Lee et al., 2003), and the prediction of fold-change robustness (Goentoro and Kirschner, 2009).

- The Tgfβ model reveals the roles of Smad nucleocytoplasmic shuttling in transmitting signaling dynamics (Schmierer et al., 2008).

We now include this information in the Introduction, to provide the readers with more evidence of the predictive power of the models.

Although the models capture systems-level behaviors, they do not necessarily reproduce fine-tuned aspects, e.g., response magnitude or the time to reach steady state, which may quantitatively vary across specific contexts. The prediction in this present study falls within this type of qualitative prediction: the models predict linear behavior, but do not necessarily capture fine-tuned aspects such as the response gain (i.e., slope of the input-output relationship).

Consistently, our aim in Figure 3 was to assess whether the measured input-output relationship was linear or not. We realize now looking at the previous version of Figure 3, that it may have come across we were making quantitative predictions, e.g., fitting the concentrations of EGF. We would like to clarify that we did not fit protein levels quantitatively, and nor does our conclusion depend on such a quantitative fit. The grey lines in Figure 3 are least squares regression lines. We have now made this more obvious in Figure 3 and the legend.

For Figure 3D specifically, we did use the model to fit the data, in the following way. Since a linear fit does not explain the data well, we tested whether the data are better captured by a nonlinear fit. As a nonlinear fit, rather than testing ad-hoc nonlinear functions, we used the full ERK model. We first fitted the gain of the model to the data (i.e., the y-range), and afterward, the only parameter varied was the strength of the negative feedback (k25) – since we mutated the negative feedback in this experiment. Using the Akaike Information Criterion (AICc), we verified that a nonlinear fit using the ERK model outperforms a least squares linear fit.

Reviewer #3:

This manuscript aims to test whether signaling pathways have linear input-output response characteristics. A simplification of existing mathematical models indicates this being the case for the canonical Wnt, ERK and TGF-β pathways, despite striking differences in core network architectures among these pathways. It is interesting to see how different complex pathways can be approximated by linear systems. The manuscript concludes by experiments (Western Blots) testing the theoretical results, as well as the breakdown of linearity as some key features of the pathways are altered.

In summary the major point here seems to be that in many of these diverse, complicated, pathways, a linear dose-response can be observed as an approximation, and this might be borne out in experimental conditions as well. In fact, not only can the linear dose response be observed, it appears rather robust to perturbations until you exceed some critical value. While these results are potentially interesting, there are major weaknesses that should be addressed before the manuscript could be considered for publication in eLife.

1) Besides the three pathways chosen, there are other pathways that are similarly well studied, both experimentally as well as mathematically. The NF-kappaB pathway is a good example. Why was the NF-kappaB pathway or any other pathway not considered? You could make a statement that proving the generality of this finding will require investigating additional pathways in the future.

We did not include the NF-κB pathway because it is analytically intractable in our hands. However, we agree with reviewer 3 that discussion of the NF-κB pathway could help motivate testing linearity in other signaling pathways, even when analytics is not possible. We numerically solved the NF-κB model using parameters that have been measured or fitted in cells. We found that, despite having strong nonlinearities e.g., sequestration, the NF-κB model shows a linear input-output relationship over a physiologically relevant input range. This further strengthens the finding of convergent linearity across multiple complex signaling pathways. We now include this finding in Discussion section and Figure 3—figure supplement 7.

Additional paragraph in Discussion section:

“It would be interesting to further probe the generality of linear signal transmission. … Besides the systems analyzed here, NF-κB is another signaling pathway that has been modeled rigorously (5, 50-51). Numerical simulations of a well-established NF-κB model (50) over the range of nuclear NF-κB translocation observed in human epithelial cells (51) reveal that the peak of the nuclear NF-κB pulse correlates linearly with ligand concentration (Figure 2—figure supplement 7).”

2) A linear input-output relationship due to negative feedback has already been described for a MAP Kinase pathway, see PMID:19079053. How are the current results novel compared to these earlier findings (besides studying a human MAPK pathway)?

To begin with, these are two distinct MAPK pathways: Fus3(MAPK) pathway in fungal yeast is coupled to G protein-coupled receptor (GPCR), whereas we are considering ERK(MAPK) pathway in mammals coupled to receptor tyrosine kinases. Even if we simply consider negative feedback broadly, the present results are novel in two ways:

1) The negative feedbacks between the fungal and animal MAPK pathways, as currently known, are convergent. Yu et al. 2008 (Richard et al., 2008) identified negative feedback by Fus3 acting on Sst2, a yeast GAP protein. It is known that ERK feedbacks on GRK, a kinase upstream from GAP, but no feedbacks have been identified on mammalian GAPs. Even the newly identified feedback actions of Fus3, on Ste5 and Ste18 (Choudhury, Baradaran-Mashinchi and Torres, 2018), are not conserved in mammals. Finally, with regard to the specific feedback we are considering here, where ERK hyper-phosphorylates Raf, there is no known homolog of Raf in yeasts.

2) The linearity we identified in the EGF/ERK model here not only requires negative feedback, but also ultrasensitivity in the Raf-MEK-ERK cascade.

These differences further strengthen the finding of convergent linearity across independently evolving signaling pathways, as well as homologous pathways that diverged 1.5 billion years ago. We now include this in Discussion section:

“Finally, linearity extends beyond metazoan signaling pathways. In the yeast pheromone sensing pathway, a homolog of the ERK cascade, transcriptional output correlates linearly with receptor occupancy (Richard et al., 2008). The linearity is mediated by negative feedback by Fus3 acting on Sst2, a feedback that is not conserved in the mammalian ERK system. These further argue for linear signal transmission as a convergent property across independently evolving signaling pathways, as well as between conserved pathways that diverged 1.5 billion years ago.”

3) The claims of linearity in this manuscript are mainly based on visual examination. However, there are rigorous ways to measure linearity, based on the L1-norm, as described in PMID:19279212 and PMID:23385595, which should be cited. It would be necessary to apply the L1-norm throughout all figures of the manuscript to make computational and experimental claims of linearity quantitative.

We apologize that some information was missing from Figure 3 that gives the impression that we assessed the linearity only visually. We have now noted more clearly in Figure 3 that we applied for each dataset least squares regression analysis (L2-norm). As there is complementarity between L1- and L2-norm, we also follow reviewer 3’s suggestion and added L1-norm analysis for all linearity assessment in the manuscript. The conclusions hold with both analyses.

(We now address Point 7 and 8 first, since these points are relevant to many other questions.)

7) Considering the importance of fold-change detection in biology, it would be important to assess the number of decades over which linearity holds for each pathway. In fact, any nonlinear function can be approximated by a linear relationship except near extremum points (this is the essence of the Taylor approximation). Are these observations more than Taylor expansion? Maybe we can tell based on the decades over which linearity holds.

As reviewer 1 also suggests, we see now that the phrase “physiological range” creates a misleading impression that we have information about the entire physiological range of the parameters across tissues and organisms. We have now described more carefully in the manuscript that we considered a physiological relevant parameter regime, i.e., the values of parameters that have been measured in some biological systems, that we specify in the main text. As linearity emerges with these measured parameters, it argues for the physiological relevance of the predicted linearity. Indeed, we observed linear input-output relationship in two human cell lines.

With this revision, it should be clear now that we do not suggest that linearity occurs in all possible physiological ranges, but in some physiological contexts. As also raised by reviewer 1 (Point 2), the ERK pathway produces ultrasensitive response in some biological contexts, as well as graded response in some other contexts. Both responses can be captured by modulating the negative feedback in the ERK model. Finally, we now include a section in the supplement entitled “Toy Model of the ERK pathway” that illustrates how ultrasensitivity can produce linearity, to further clarify their seemingly contradictory relationship.

8) The data in Figure 3 should be expanded into the domains of nonlinearity/saturation and then the range of linearity should be indicated on top of such curves (linearity measured using the L1-norm). This should be done for all panels. In fact, without showing the full curve (including the nonlinear parts), panel 3C does not indicate that linearity is lost or that its range shrinks.

As also suggested by reviewer 2, we now explain more explicitly that the predicted linearity extends over the dynamic range of the model, theoretical or observed (Figure 2—figure supplement 2).

Across these three natural pathways, the dynamic range may not appear as impressive as e.g., some synthetic circuits that can operate across multiple orders of magnitudes. And yet, this is one reason, we think, fold-change detection may be important, as it allows the system to continually rescale its dynamic range to present stimulation level, and maintains sensitivity to respond (as discussed in ref. (Olsman and Goentoro, 2016)).

4) For deriving formula [4], apha/u>>1 is necessary. However, based on the SI, α/u=11, and it is actually smaller for some of the range plotted in Figure 2. Then the linearity arises from 'α / u >> 1 + γ'. One could/should explore when α = (1+γ)/u * S where S is a scaling factor set to say, 0.9, 1, 1.1, 2, 10, and 100. At 10 and 100, one might expect the condition to hold, but at S around 1 the simplification probably fails. Do we then see lack of linearity? One can then re-arrange the terms for each of the other terms (u or γ) and explore similarly. The same argument can be made for the Erk and Tgf[β] pathways.

We thank reviewer 3, we now include this analysis as Figure 2—figure supplement 5.

5) How robust is linearity? To answer this question, two modifications are needed: (i) extend the range of u for each plot in Figure 2, showing the nonlinearity of the curves, then indicate the linear range on such nonlinear curves, e.g. with a different color; (ii) Change a parameter and plot a family of curves, indicating the linear range (determined based on the L1-norm) on each individual curve. This can be done for a couple of parameters, but it is especially important for the parameters mediating interactions that cause linearity to break down (described in subsection “Linearity in the Wnt and ERK pathways is modulated by perturbation to parameters”).

We agree this are a good analysis to include. For the input range, as discussed in Point 8 above, the linearity extends the dynamic range of the models. For parameter variation, we now include this analysis as a supplement figure. In these simulations, we plot input-output relationship over the entire dynamic range of the models in unperturbed system (as discussed in Point 8). Shaded in purple is the region of the response with L1-norm < 0.1. The analysis shows that linearity occurs through a considerable range of parameter values.

6) If linear responses indeed help cellular signal processing as a result of convergent evolution then linearity should probably occur in single cells that process inputs. Unfortunately, the Western Blots in Figure 3 only test linearity at the population average level. The linearity of the average does not imply linearity in single cells. So, are single cell input-output response characteristics linear? This should be checked at least for one signaling pathway, possibly by new experiments or otherwise using data by others, see for example PMID:25504722. In fact, recent papers indicate that some pathways' responses are dynamic/oscillatory, noisy or bimodal at the single cell level. All of this should be addressed/discussed. Where is then the linearity?

We agree with the reviewer that discussing this point is necessary. While quantitative assessment of single-cell dynamics is now increasingly done, double quantitation of input-output in single cells is still technically challenging. It is even more challenging in our study because the inputs and outputs here are not simply protein levels, but include phosphorylated proteins (e.g., pLRP, dpERK). On the advice of the reviewer, we examined existing single-cell data from ref. (Selimkhanov et al., 2014), kindly provided by Roy Wollman, which utilizes a FRET sensor to measure ERK activity. The dataset offers useful insights into the dynamics of single-cell ERK response (as discussed in their paper); however, the dynamic range of the sensor was not large enough to distinguish between dose-response models (<2-fold, compared to 12-fold for our Western Blot data). Another method, quantitative dual-immunofluorescence, requires having at least two good linear antibodies against each target of interest to calibrate the linearity of the antibodies (e.g., what was done in ref. (Cheong et al., 2011)), and unfortunately, which are not available for the targets of interest here.

Despite these technical limitations, we can still address reviewer’s concern in 2 ways. First, we can qualitatively assess existing data from single cells for the plausibility of linearity in some biological contexts. Linear behavior means that single cells responds to ligand in a graded manner. Even though there are reports of oscillatory or bimodal response in signaling pathways, there are also multiple observations where signaling pathways show graded response in single cells. We provide some of these examples for the three pathways in Supplementary file 4.

We include this information in Discussion section:

“It would be interesting to further probe the generality of linear signal transmission. Linear behavior requires that single cells responds to ligand in a graded manner. Although there are reports of oscillatory or bimodality in signaling pathways, there are also multiple observations across biological contexts of single cells responding to ligand in a graded manner (Table S4).”

Second, the ERK pathway, in particular, is sometime observed to show ultrasensitive response. We confirmed that the particular cells used in this study, non-small cell lung carcinoma, have been observed to show graded response at single-cell level (measured in live cells at the level of total ERK, (Cheong et al., 2011)), with no evidence for bimodality. We further confirmed qualitatively the graded response in these cells using immunofluorescence against dpERK (now included as Figure 3—figure supplement 4) (note as described earlier, we cannot quantify this immunofluorescent signal as we cannot calibrate if the antibody is linear in immunofluorescence assay).

9) The data in Figure 3 should be expanded into the domains of nonlinearity/saturation and then the range of linearity should be indicated on top of such curves (linearity measured using the L1-norm). This should be done for all panels. In fact, without showing the full curve (including the nonlinear parts), panel 3C does not indicate that linearity is lost or that its range shrinks.

As also requested by reviewer 2, we have now included more measurements to further confirm the dynamic range of the cells used in the experiments. The revised Figure 3 is included below, where data points beyond saturation are shown in grey, and linearity is assessed in blue data points.

10) What ranges on the theoretical plots do the experimental results correspond to? Theory and experiment should be better connected.

We have now made the theory-experiment connection clearer by making it more explicit in text and Figure 2—figure supplement 2 that the models predict linearity across the dynamic range of the systems (as detailed in Point 8). Corresponding to the model prediction, we measured the entire dynamic range of the cells and observed that linearity governs the entire range until saturation (as detailed in Point 9).

11) Experimental verification is not trivial in these systems. First, the "quantitative' Western blots are not truly quantitative, despite ensuring that fluorescent antibodies are within linear range, etc. Western blots, by their very nature, are notoriously noisy, and therefore require a lot of "normalization" which might artifactually introduce the appearance of linearity. If, for example, the denominator is very large for normalization, the signal can appear more linear. Say, if you compare y=x and y=x^2, but then divide by a large factor K such that within the range of 'x' being explored K>>y, then both equations will be nearly 0 and will likely show a "linear" response (also true according to the Taylor expansion near 0).

To show that the linearity is not a caveat of normalization in the measurements, here are some raw data without loading control normalization, in experiments where the loading control across the samples came out fairly uniform, <10%. (This can sometime happen because we processed samples in our protocol as uniformly as possible). The observed linearity holds even without normalization.

12) The variability in the Western blots is borne out in Figure 3A where the various 'Wnt' doses have very variable pLRP5/6 response (input). Likewise, in Figure 3B, for EGF. These few points do not necessarily demonstrate linearity other than the fact that the linear approximation can be applied to just about any relationship.

The variability of our Western blot measurement is within 10%, arguing for the precision of our measurements (demonstrated in figure below, now included as Figure 3—figure supplement 7).

While we can tightly control technical variability, there are inevitable biologic differences, e.g., each plate of cells may have slightly different cell confluency, cell contact, cell cycle state and therefore slightly different receptivity to ligands. In addition, measurements of ligand concentrations may also have uncertainty. All this may influence response amplitude (or gain) across experiments, which is also a highly fine-tuned parameter in the models. Despite such variability in response gain across experiments, the linear response is highly reproducible across experiments. The following figure shows some measurements from independent experiments. Moreover, where we can measure the internal proxy for ligand activation (phopho-LRP), regardless of uncertainty in externally measured Wnt dose, it is notable that phospho-LRP5/6 remains linearly correlated with β-catenin.

12) (ii) In addition, comparing Figure 3A with Figure 3D where linearity breaks down with a Raf mutation preventing feedback, it still appears linear until 3ng/mL, a dose that is not reported for the wild-type (Figure 3A). In fact, examining both plots between 0 to 2ng/mL, one could conclude instead that both systems are 'linear'.

We agree that the way we currently use the term ‘linearity’ in the text can be confusing, since any function would trivially have a linear regime for some small range. Following reviewer 2’s suggestion, we have now made it more clearly in the revised manuscript that the predicted linearity spans the dynamic range of the models. Consistently, in testing the prediction in the ERK pathway, we analyzed linearity throughout the dynamic range of the cells (from 0-4 ng/mL EGF). We now include in the main text a data point for 3 ng/mL EGF in wt cells as well. Performing L1-norm analysis over 0-4 ng/mL EGF, we find for wt cells, L1-norm = 0.03, and for mutant cells, L1-norm = 0.26. Thus, one may conclude that upon mutation in the negative feedback linearity was lost or reduced, either way confirming that the negative feedback is necessary to maintain linearity.

13) Same with the Wnt pathway. Here it is already recognized that linearity was not lost, but this is attributed to side effects of the drug used – something that could be tested perhaps with a different drug or a mutation similar to the Erk pathway analysis.

Indeed, in the Wnt pathway, we not only observed linearity, but also that it is unexpectedly more strongly maintained in the cells than the model predicts. This is not a side effect of the drug, but rather a feature of the dual, incoherent role of GSK3β in phosphorylating βcatenin for degradation and phosphorylating LRP for the pathway activation (phospho-LRP inhibits β-catenin degradation). Encoding such an incoherent feedforward loop in the model does capture the maintenance of linearity. We now include the model in Figure 2—figure supplement 4 and add this discussion in the main text (Discussion section). Thus, the experiments not only confirm linear signal transmission in the Wnt pathway, but also reveal a hitertho unknown role of the GSK3β incoherent feedforward loop in maintaining the linearity.

https://doi.org/10.7554/eLife.33617.035

Article and author information

Author details

  1. Harry Nunns

    Division of Biology and Biological Engineering, California Institute of Technology, Pasadena, United States
    Contribution
    Conceptualization, Data curation, Software, Investigation, Methodology, Writing—original draft, Writing—review and editing
    For correspondence
    hnunns@caltech.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9669-0039
  2. Lea Goentoro

    Division of Biology and Biological Engineering, California Institute of Technology, Pasadena, United States
    Contribution
    Conceptualization, Funding acquisition, Methodology, Writing—original draft, Writing—review and editing
    For correspondence
    goentoro@caltech.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3904-0195

Funding

James S. McDonnell Foundation (220020365)

  • Lea Goentoro

National Science Foundation (NSF.145863)

  • Lea Goentoro

National Institutes of Health (5T32GM007616-37)

  • Harry Nunns

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We would like to thank Rob Oania for providing advice on experiments, Michael Abrams, Christopher Frick, Kibeom Kim, and Noah Olsman for comments on the manuscript, and Michael Elowitz and Richard Murray for discussions on the study.

Senior Editor

  1. Aviv Regev, Broad Institute of MIT and Harvard, United States

Reviewing Editor

  1. Wenying Shou, Fred Hutchinson Cancer Research Center, United States

Reviewers

  1. Wenying Shou, Fred Hutchinson Cancer Research Center, United States
  2. Steven S Andrews, Fred Hutchinson Cancer Research Center, United States

Publication history

  1. Received: November 16, 2017
  2. Accepted: September 9, 2018
  3. Accepted Manuscript published: September 17, 2018 (version 1)
  4. Accepted Manuscript updated: September 18, 2018 (version 2)
  5. Accepted Manuscript updated: September 19, 2018 (version 3)
  6. Version of Record published: October 25, 2018 (version 4)

Copyright

© 2018, Nunns et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,610
    Page views
  • 408
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Cell Biology
    Eutteum Jeong et al.
    Research Article Updated
    1. Cell Biology
    Sam A Menzies et al.
    Research Communication