Signaling pathways as linear transmitters

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record published: October 25, 2018 (This version)
Accepted Manuscript updated: September 19, 2018 (Go to version)
Accepted Manuscript updated: September 18, 2018 (Go to version)
Accepted Manuscript published: September 17, 2018 (Go to version)
Accepted: September 9, 2018
Received: November 16, 2017

1. Related to
Signaling Systems: Transferring information without distortion

Steven S Andrews, Roger Brent, Gábor Balázsi

Insight Oct 25, 2018
Further reading

Abstract
Introduction
Results
Discussion
Materials and methods
Appendix 1
Data availability
References
Decision letter
Author response
Article and author information
Metrics

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.

Figure 1

Download asset Open asset

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):

{[β c a t]}_{s s} = K_{17} \cdot \frac{1 - γ + \frac{α}{u}}{2} (\sqrt{1 + \frac{4 γ}{(1 - γ + {\frac{α}{u}}^{2})}} - 1)

α = \frac{k_{4} k_{6} k_{9} v_{14} \cdot G S K 3_{t o t} \cdot A P C_{t o t}}{k_{5} k__{6} K_{7} K_{8} k_{13} k_{15}}

γ = \frac{v_{12}}{k_{13} K_{17}}

Figure 2 with 7 supplements see all

Download asset Open asset

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 k₁₀ and k₁₁ 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 (W n t)$ is the rate by which Dishevelled/Dvl inhibits the destruction complex upon Wnt ligand activation, where $k_{3}$ and $k_{- 6}$ 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

Figure 2—source code 1: https://doi.org/10.7554/eLife.33617.011
Download elife-33617-fig2-code1-v4.zip

where the input function $u = u (W n t)$ is the rate of inhibition of the destruction complex (DC) via Dishevelled/Dvl, a function of ligand-receptor activation. As illustrated in Figure 2A, K_i’s are equilibrium dissociation constants, k_i’s are rate constants, and v_i’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, $α / u ≫ 1$ 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 $α / u ≫ 1 + γ$ , Equation 1 simplifies to

{[β c a t]}_{s s} \approx K_{17} \frac{γ}{α} 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 ( $u \sim 6.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 $α / u ≫ 1 + γ$ (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 ${[d p E R K]}_{s s}$ and the input to the phosphorylation cascade $u$ :

{[d p E R K]}_{s s} = \frac{α}{β} \cdot (\frac{R a f_{t o t}}{[p R a f]_{s s}}) - 1 - \frac{γ}{α} \cdot u - \frac{δ}{β}

α = \frac{k_{3} \cdot (k_{8} + k_{b 7})}{k_{7} \cdot [P 1]_{s s} \cdot k_{8}} + \cdot \cdot \cdot

β = \frac{k_{25} \cdot (k_{30} + k_{b 29} + k_{29} \cdot [P 4]_{s s})}{k_{29} \cdot [P 4]_{s s} \cdot k_{30}} + \cdot \cdot \cdot

γ = \frac{k_{3} \cdot (k_{8} + k_{b 7}) \cdot k_{9} \cdot [M E K]_{s s}}{k_{7} \cdot [P 1]_{s s} \cdot k_{8} \cdot k_{10}} + \cdot \cdot \cdot

δ = \frac{k_{26} + k_{b 25}}{k_{26}} + \dots

Detailed derivations of Equation 5 are presented in Appendix 1. The input $u = u (E G F)$ 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, A29–A31). The dimensionless group $α \cdot u$ relates to the amount of free, phosphorylated Raf ( $α$ , blue-shaded in Figure 2B), $β \cdot {[d p E R K]}_{s s}$ 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 $γ \cdot 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 ${[p R a f]}_{s s}$ , 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 $β \cdot {[d p E R K]}_{s s} ~ (α + γ) \cdot 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 ${[d p E R K]}_{s s}$ with respect to ${[p R a f]}_{s s}$ , ultrasensitivity entails that $K ≫ 1$ . 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 $β {\cdot [d p E R K]}_{s s} ~ (α + γ) \cdot u + δ$ and $K ≫ 1$ , the negative feedback holds the level of pRaf constant ( ${[p R a f]}_{s s} \approx R_{s}$ , details in Appendix 1). With these two features, strong negative feedback and ultrasensitivity, dpERK becomes a linear function of the input $u$ :

{[d p E R K]}_{s s} \approx \frac{α}{β} \cdot \frac{R a f_{t o t}}{R_{s}} \cdot u - \frac{δ}{β}

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:

{[S 24 n]}_{s s} = a \cdot \frac{α \cdot u}{(α + γ) \cdot u + β} S 2_{t o t}

α = \frac{a \cdot (k_{o n} [S 4 n]_{s s} + a \cdot k_{e x 2})}{k_{o f f}} + \cdot \cdot \cdot

β = \frac{P P a s e \cdot k_{d e p h o s}}{k_{p h o s} \cdot R_{t o t} \cdot \frac{k_{e x 2}}{a \cdot k_{e x 2} + k_{i n 2}}} + \cdot \cdot \cdot

γ = a \cdot (a \cdot k_{e x 2} + P P a s e \cdot k_{d e p h o s}) (\frac{1}{a \cdot k_{e x 2}} + \frac{1}{C I F \cdot k_{i n 2}}) + \cdot \cdot \cdot

In Equation 11, the input function $u = u (T g f β)$ 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): $α \cdot 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 $γ \cdot 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 ${[S 4 n]}_{s s}$ 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 ( $u \sim 0.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$ , $α \cdot u ~ 1.5$ , and $γ \cdot 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 ( $β ≫ (α + γ) \cdot 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 $β ≫ (α + γ) \cdot u$ , the first term in the denominator of Equation 11 is small, and concentration of nuclear Smad complex becomes a linear function of input:

{[S 24 n]}_{s s} \approx a \cdot \frac{α \cdot S 2_{t o t}}{β} \cdot 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 $β ~ (α + γ) \cdot 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

Download asset Open asset

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 r² is the coefficient of correlation where r² = 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 (k₂₅) 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

Figure 3—source data 1: https://doi.org/10.7554/eLife.33617.022
Download elife-33617-fig3-data1-v4.zip

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. $α / u ≫ 1 + γ$ ). 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).

Figure 4

Download asset Open asset

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 = a \cdot u$ , 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 $F C = Y_{s t i m u l a t e d} / Y_{b a s a l}$ 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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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, $r^{2}$ .

Akaike information criterion

Request a detailed protocol

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 (A I C)$ 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

Download asset Open asset

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

Download asset Open asset

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,

\dot{[A]} = - k_{1} [A] \cdot [D]^{2} + k_{3} [C] = 0

\dot{[B]} = k_{1} [A] \cdot [D]^{2} - k_{2} [B] = 0

\dot{[C]} = k_{2} [B] - k_{3} [C] = 0

\dot{[D]} = - 2 k_{1} [A] \cdot D^{2} + k_{2} [B] + k_{3} [C] = 0

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

T_{1} = [A] + [B] + [C]

T_{2} = 2 [B] + [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:

(\begin{matrix} - k_{1} {[D]}_{s s}^{2} & 0 & k_{3} \\ k_{1} {[D]}_{s s}^{2} & - k_{2} & 0 \\ 0 & k_{2} & - k_{3} \end{matrix}) (\begin{matrix} {[A]}_{s s} \\ {[B]}_{s s} \\ [{C]}_{s s} \end{matrix}) = 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:

{[A]}_{s s} = c \cdot k_{2} k_{3}

{[B]}_{s s} = c \cdot k_{1} k_{3} D_{s s}^{2}

{[C]}_{s s} = c \cdot k_{1} k_{2} D_{s s}^{2}

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]}_{s s}$ . For instance, the solution for ${[C]}_{s s}$ is below.

{[C]}_{s s} = \frac{k_{1} k_{2} {[D]}_{s s}^{2}}{k_{2} k_{3} + k_{1} {[D]}_{s s}^{2} (k_{2} + k_{3})} T_{1}

The solutions for ${[A]}_{s s}$ , ${[B]}_{s s}$ , and ${[C]}_{s s}$ derived from the variable elimination technique still depend on ${[D]}_{s s}$ . If we plug in the solutions for the cut species, we can obtain polynomial equations for the remaining species (in this case ${[D]}_{s s}$ ), 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, $k_{2} k_{3}$ , $k_{1} k_{3} {[D]}_{s s}^{2}$ , and $k_{1} k_{2} {[D]}_{s s}^{2}$ represent the un-normalized fraction of $T_{1}$ that exists as A, B, and C, respectively. The normalization factor for these fractions is $c / T_{1}$ , 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 $T_{1}$ . For instance, increasing the value of $k_{1}$ will increase the amount of $T_{1}$ that exists as $B$ and $C$ , while necessarily decreasing the amount of $A$ (assuming ${[D]}_{s s}$ 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):

[β {c a t]}_{s s} = K_{17} \frac{1 - γ + \frac{α}{u (W n t)}}{2} (\sqrt{1 + \frac{4 γ}{{(1 - γ + \frac{α}{u (W n t)})}^{2}}} - 1)

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

α = \frac{k_{4} \cdot k_{6} \cdot k_{9} \cdot v_{14} \cdot G S K 3_{t o t} \cdot A P C_{t o t}}{k_{5} \cdot k__{6} \cdot K_{7} \cdot K_{8} \cdot k_{13} \cdot k_{15}}

γ = \frac{v_{12}}{k_{13} \cdot K_{17}}

u (W n t) = 1 + \frac{k_{3} \cdot D v l_{t o t}}{k__{6}} \cdot \frac{k_{1} \cdot W n t}{k_{2} + k_{1} \cdot W n t}

The input function $u = u (W n t)$ 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

u (W n t = 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:

\sqrt{1 + \in} \approx 1 + \frac{1}{2} \in, \in ≪ 1

\in= \frac{4 γ}{{(1 - γ + \frac{α}{u})}^{2}}

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

[β {c a t]}_{s s} \approx K_{17} \frac{γ}{α} 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. $k_{25}$ characterizes the negative feedback from dpERK to unphosphorylated Raf, and $k_{27}$ 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 $k_{25}$ , and set $k_{27}$ 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,

{[d p E R K]}_{s s} = g ({[p R a f]}_{s s})

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

{[p R a f]}_{s s} = f (u, {[d p E R K]}_{s s})

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:

{R a f, R a f : R a s G T P, p R a f, p R a f : P 1, M E K : p R a f, p M E K : p R a f, R a f : p p E R K, R a f i, R a f i : P 4}

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:

{P 1, M E K, p M E K, d p E R K, P 4}

With this, we derive the expression for ${[p R a f]}_{s s}$ ,

{[p R a f]}_{s s} = \frac{α \cdot u}{β \cdot [d p E R K]_{s s} + (α + γ) \cdot u + δ} \cdot R a f_{t o t}

where the parameter groups are:

α = \frac{k_{3} \cdot (k_{8} + k_{b 7})}{k_{7} \cdot [P 1]_{s s} \cdot k_{8}} + \cdot \cdot \cdot

β = \frac{k_{25} \cdot (k_{30} + k_{b 29} + k_{29} \cdot [P 4]_{s s})}{k_{29} \cdot [P 4]_{s s} \cdot k_{30}} + \cdot \cdot \cdot

γ = \frac{k_{3} \cdot (k_{8} + k_{b 7}) \cdot (k_{9} \cdot [M E K]_{s s})}{k_{7} \cdot [P 1]_{s s} \cdot k_{8} \cdot k_{10}} + \cdot \cdot \cdot

δ = \frac{k_{26} + k_{b 25}}{k_{26}} + \dots

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

\begin{array}{ll} α = & (k_{3} \cdot (k_{8} + k_{b 7}) \cdot (k_{10} + k_{b 9}) \cdot (k_{12} + k_{b 11}) \cdot (k_{26} + k_{b 25})) / \\ ({[P 1]}_{s s} \cdot k_{7} \cdot k_{8} \cdot k_{10} \cdot k_{12} \cdot k_{26}) \end{array}

\begin{array}{ll} β = & (k_{25} \cdot (k_{4} + k_{b 3}) \cdot (k_{10} + k_{b 9}) \cdot (k_{12} + k_{b 11}) \cdot (k_{26} \cdot k_{30} + k_{26} \cdot k_{b 29} + {[P 4]}_{s s} \cdot k_{26} \\ \cdot k_{29} + {[P 4]}_{s s} \cdot k_{29} \cdot k_{30})) / ({[P 4]}_{s s} \cdot k_{4} \cdot k_{10} \cdot k_{12} \cdot k_{26} \cdot k_{29} \cdot k_{30}) \end{array}

\begin{array}{ll} γ = & (k_{3} \cdot (k_{26} + k_{b 25}) \cdot (k_{4} \cdot k_{8} \cdot k_{10} \cdot k_{12} + k_{4} \cdot k_{8} \cdot k_{10} \cdot k_{b 11} + k_{4} \cdot k_{8} \cdot k_{12} \cdot k_{b 9} + k_{4} \\ \cdot k_{10} \cdot k_{12} \cdot k_{b 7} + k_{4} \cdot k_{8} \cdot k_{b 9} \cdot k_{b 11} + k_{4} \cdot k_{10} \cdot k_{b 7} \cdot k_{b 11} + k_{4} \cdot k_{12} \cdot k_{b 7} \\ \cdot k_{b 9} + k_{4} \cdot k_{b 7} \cdot k_{b 9} \cdot k_{b 11} + {[M E K]}_{s s} \cdot k_{4} \cdot k_{8} \cdot k_{9} \cdot k_{12} + {[M E K]}_{s s} \cdot k_{4} \cdot k_{8} \\ \cdot k_{9} \cdot k_{b 11} + {[M E K]}_{s s} \cdot k_{4} \cdot k_{9} \cdot k_{12} \cdot k_{b 7} + {[M E K]}_{s s} \cdot k_{4} \cdot k_{9} \cdot k_{b 7} \cdot k_{b 11} \\ + {[P 1]}_{s s} \cdot k_{4} \cdot k_{7} \cdot k_{10} \cdot k_{12} + {[P 1]}_{s s} \cdot k_{7} \cdot k_{8} \cdot k_{10} \cdot k_{12} + {[P 1]}_{s s} \cdot k_{4} \cdot k_{7} \cdot k_{10} \\ \cdot k_{b 11} + {[P 1]}_{s s} \cdot k_{4} \cdot k_{7} \cdot k_{12} \cdot k_{b 9} + {[P 1]}_{s s} \cdot k_{7} \cdot k_{8} \cdot k_{10} \cdot k_{b 11} + {[P 1]}_{s s} \cdot k_{7} \\ \cdot k_{8} \cdot k_{12} \cdot k_{b 9} + {[P 1]}_{s s} \cdot k_{4} * k_{7} \cdot k_{b 9} \cdot k_{b 11} + {[P 1]}_{s s} \cdot k_{7} \cdot k_{8} \cdot k_{b 9} \cdot k_{b 11} \\ + k_{4} \cdot k_{8} \cdot k_{10} \cdot k_{11} \cdot {[p M E K]}_{s s} + k_{4} \cdot k_{8} \cdot k_{11} \cdot k_{b 9} \cdot {[p M E K]}_{s s} + k_{4} \cdot k_{10} \\ \cdot k_{11} \cdot k_{b 7} \cdot {[p M E K]}_{s s} + k_{4} \cdot k_{11} \cdot k_{b 7} \cdot k_{b 9} \cdot {[p M E K]}_{s s})) / \\ ({[P 1]}_{s s} \cdot k_{4} \cdot k_{7} \cdot k_{8} \cdot k_{10} \cdot k_{12} \cdot k_{26}) \end{array}

δ = ((k_{4} + k_{b 3}) \cdot (k_{10} + k_{b 9}) \cdot (k_{12} + k_{b 11}) \cdot (k_{26} + k_{b 25})) / (k_{4} \cdot k_{10} \cdot k_{12} \cdot k_{26})

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, $α \cdot u$ relates to the amount of free, phosphorylated Raf since $α \cdot u / ((α + γ) \cdot u + β {[d p E R K]}_{s s} + δ)$ is the fraction of Raf present as pRaf. Thus, as $α \cdot u$ increases relative to $γ \cdot u + β {[d p E R K]}_{s s} + δ$ , the amount of pRaf also increases.

We can define three subpopulations of Raf: Raf inhibited by dpERK, $[R_{i}]$ ; Raf activated by RasGTP (input), $[R_{a}]$ ; and unphosphorylated Raf $[R_{n}]$ . Specifically:

[R_{i}] = [R a f : d p E R K] + [R a f i] + [R a f : P 4]

[R_{a}] = [p R a f] + [p R a f : P 1] + [M E K : p R a f] + [p M E K : p R a f] + R a f : R a s G T P]

[R_{n}] = [R a f]

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

{[R_{i}]}_{s s} = \frac{β \cdot [d p E R K]_{s s}}{((α + γ) \cdot u + β [d p E R K]_{s s} + δ)} R a f_{t o t}

{[R_{a}]}_{s s} = \frac{γ \cdot u}{((α + γ) \cdot u + β [d p E R K]_{s s} + δ)} R a f_{t o t} + [p R a f]_{s s}

{[R_{n}]}_{s s} = \frac{δ}{((α + γ) \cdot u + β [d p E R K]_{s s} + δ)} \cdot R a f_{t o t}

Thus, in the same sense that $α \cdot u$ relates to the amount of free phosphorylated Raf, $β \cdot {[d p E R K]}_{s s}$ relates to the amount of inhibited Raf, $γ \cdot 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 ${[d p E R K]}_{s s} = g ({[p R a f]}_{s s})$ 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 $d p E R K$ (from 10% to 90% of max, Figure 2—figure supplement 1B-C).

We therefore approximate ${[p R a f]}_{s s}$ by a value $R_{s}$ 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 ${[d p E R K]}_{s s}$ becomes a linear function of input:

{[d p E R K]}_{s s} \approx \frac{α}{β} \cdot (\frac{R a f_{t o t}}{R_{s}} - 1 - \frac{γ}{α}) \cdot u - \frac{δ}{β}

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

\frac{α}{β} \cdot \frac{R a f_{t o t}}{R_{s}} = 140

\frac{α}{β} (1 + \frac{γ}{α}) = 13

We can neglect the second term, yielding:

{[d p E R K]}_{s s} \approx \frac{α_{1}}{β} \cdot \frac{R a f_{t o t}}{R_{s}} \cdot u - \frac{δ}{β}

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

{[d p E R K]}_{s s} = g ({[p R a f]}_{s s})

{[{p R a f}_{t o t}]}_{s s} = f (u, {[d p E R K]}_{s s})

We can now derive a general expression for the relative change of ${[p R a f]}_{s s}$ with respect to a relative change in $u$ . We use the notation $d \hat{x} = d \ln x = d x / x$ .

\frac{d \hat{f}}{d \hat{u}} = \frac{\partial f}{\partial u} \cdot \frac{u}{f} \cdot {(1 - \frac{\partial f}{\partial [d p E R K]_{s s}} \cdot \frac{d g}{d [p R a f]_{s s}})}^{- 1}

Next, we define the response coefficient $K$ between ${[d p E R K]}_{s s}$ and ${[p R a f]}_{s s}$ :

K ≜ \frac{d g}{d [p R a f]_{s s}} \cdot \frac{[p R a f]_{s s}}{[d p R a f]_{s s}}

From Equation A21, we get the partial derivatives:

\frac{\partial f}{\partial u} = \frac{f}{u} \cdot \frac{β [d p E R K]_{s s} + δ}{β [d p E R K]_{s s} + (α + γ) u + δ}

\frac{\partial f}{\partial {[d p E R K]}_{s s}} = - f \cdot \frac{β}{β [d p E R K]_{s s} + (α + γ) u + δ}

Using these two equations, we find that:

\frac{d \hat{f}}{d \hat{u}} = \frac{β {[d p E R K]}_{s s} + δ}{(1 + k) \cdot β [d p E R K]_{s s} + α u + δ}

When $K ≫ 1$ and $β {[d p E R K]}_{s s} ~ (α + γ) \cdot u + δ$ , we see that

\frac{d \hat{f}}{d \hat{u}} \approx K^{- 1}

Therefore, ${[p R a f]}_{s s}$ 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 ${[d p E R K]}_{s s}$ becomes a linear function of input.

\frac{d g \hat{-} g_{0}}{d \hat{u}} \approx 1; g_{0} = - \frac{δ}{β}

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:

An input $u$ increases the amount of species $R$ , which in turn influences $E$ as $E = g (R)$ . There is negative feedback from $E$ to $R$ , which in the limit of strong negative feedback is inversely proportional to $E$ .
Next, we specify the function $g (R)$ such that $K = K_{0}$ , where $K$ is the relative change of $E$ with respect to $R$ . As $K_{0}$ increases, therefore, the function $g (R)$ becomes more ultrasensitive.

Solving for $E$ , we see that in the limit of $K = K_{0} ≫ 1$ , $E$ becomes a linear function of $u$ , and $R$ is held constant at R_s.

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 $K ≫ 1$ . 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

Download asset Open asset

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:

{S 2 c, p S 2 c, S 24 c, S 2 n, p S 2 n, S 24 n}

which is subject to the conservation equation:

([S 2 c] + [p S 2 c] + [S 24 c]) + \frac{1}{a} ([S 2 n] + [p S 2 n] + [S 24 n]) = S 2_{t o t}

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

{S 4 c, S 4 n}

Using this relationship, we derive an expression for the nuclear Smad complex $(S 24 n$ ) at steady-state,

{[S 24 n]}_{s s} = \frac{a \cdot α \cdot u}{(α + γ) \cdot u + β} S 2_{t o t}

where the parameter groups are:

α = \frac{a \cdot (k_{o n} [S 4 n]_{s s} + a \cdot k_{e x 2})}{k_{o f f}} + \cdot \cdot \cdot

β = P P a s e \cdot \frac{k_{d e p h o s}}{k_{p h o s} \cdot R_{t o t} \cdot \frac{k_{e x 2}}{a \cdot k_{e} x 2} + k_{i n 2}} + \cdot \cdot \cdot

γ = a \cdot (a \cdot k_{e x 2} + P P a s e \cdot k_{d e p h o s}) (\frac{1}{a \cdot k_{e x 2}} + \frac{1}{C I F \cdot k_{i n 2}}) + \cdot \cdot \cdot

Here the input function $u = u (T g f β)$ 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 ${[S 4 c]}_{s s}$ and ${[S 4 n]}_{s s}$ calculated for $u = 0$ ). All calculations for the paper use these truncated parameter groups. The complete parameter groups are written below:

\begin{array}{ll} α = & (a \cdot ({[S 4 n]}_{s s} \cdot k_{o f f} + C I F \cdot {[S 4 n]}_{s s} \cdot k_{i n 2} + C I F \cdot P P a s e \cdot {[S 4 c]}_{s s} \cdot k_{d e p h o s} + C I F \\ \cdot {[S 4 c]}_{s s} \cdot {[S 4 n]}_{s s} \cdot k_{o n} + C I F \cdot {[S 4 c]}_{s s} \cdot a \cdot k_{e x 2})) / (C I F \cdot {[S 4 c]}_{s s} \cdot k_{o f f}) \end{array}

\begin{array}{ll} β = & (P P a s e \cdot k_{d e p h o s} \cdot (k_{i n 2} + a \cdot k_{e x 2}) \cdot (k_{o f f} + C I F \cdot k_{i n 2} + C I F \cdot {[S 4 c]}_{s s} \cdot k_{o n})) / \\ (C I F \cdot R_{t o t} \cdot {[S 4 c]}_{s s} \cdot k_{e x 2} \cdot k_{o n} \cdot k_{p h o s}) \end{array}

\begin{array}{ll} γ = & ((P P a s e \cdot k_{d e p h o s} + a \cdot k_{e x 2}) \cdot (k_{i n 2} \cdot k_{o f f} + C I F \cdot k_{i n 2}^{2} + a \cdot k_{e x 2} \cdot k_{o f f} + C I F \cdot {[S 4 c]}_{s s} \\ \cdot k_{i n 2} \cdot k_{o n} + C I F \cdot a \cdot k_{e x 2} \cdot k_{i n 2} + {[S 4 c]}_{s s} \cdot a \cdot k_{e x 2} \cdot k_{o n})) / \\ (C I F \cdot {[S 4 c]}_{s s} \cdot k_{e x 2} \cdot k_{i n 2} \cdot k_{o n}) \end{array}

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, $α \cdot u$ relates to the amount of nuclear Smad complex since $α \cdot u / ((α + γ) \cdot u + β)$ is the fraction of Smad2 present as S24n. Thus, as $α \cdot u$ increases relative to $γ \cdot u + β$ , the amount of S24n also increases.

By definition, the parameter groups $β$ and $γ \cdot 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:

{[S 2 c]}_{s s} + \frac{1}{a} {[S 2 n]}_{s s} = \frac{β + δ \cdot u}{β + (α + γ) \cdot u} S 2_{t o t}

δ = P P a s e \cdot k_{d e p h o s} \cdot \frac{k_{o f f} + C I F \cdot k_{i n 2} + C I F \cdot [S 4 c]_{s s} \cdot k_{o n}}{C I F \cdot [S 4 c]_{s s} \cdot k_{e x 2} \cdot k_{o n}}

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

{[S 2 c]}_{s s} + \frac{1}{a} {[S 2 n]}_{s s} \approx \frac{β}{β + (α + γ) \cdot u} S 2_{t o t}

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

{[p S 2 c]}_{s s} + {[S 24 c]}_{s s} + \frac{1}{a} {[p S 2 n]}_{s s} = \frac{(γ - δ) \cdot u}{((α + γ) \cdot u + β)} \cdot S 2_{t o t}

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

α \cdot u = 3.1

γ \cdot u = 1.3

β = 46

where we have used a non-saturating input ( $u = 0.2$ ). Therefore, with measured parameters, $β ≫ (α + γ) \cdot u$ . With this, the denominator in the $[S 4 {n]}_{s s}$ equation simplifies, and the concentration of Smad complex becomes a linear function of the input:

{[S 24 n]}_{s s} \approx \frac{α \cdot S 2_{t o t}}{β} \cdot u

Appendix 1—table 1

Parameters, variables, and equations of the Wnt model.

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

Parameter	Label	Value
Activation rate of Disheveled/Dvl by Wnt	$k_{1}$	$0.182$	$m i n^{- 1}$
Inactivation rate of Dvl	$k_{2}$	$1.82 \cdot 10^{- 2}$	$m i n^{- 1}$
Dissociation of destruction complex (DC) by active Dvl	$k_{3}$	$5.00 \cdot 10^{- 2}$	${n M}^{- 1} m i n^{- 1}$
Phosphorylation of DC	$k_{4}$	$0.267$	$m i n^{- 1}$
Dephosphorylation of DC	$k_{5}$	$0.133$	$m i n^{- 1}$
Forward rate for DC binding	$k_{6}$	$9.09 \cdot 10^{- 2}$	$n M^{- 1} m i n^{- 1}$
Reverse rate for DC binding	$k_{- 6}$	$0.909$	$m i n^{- 1}$
Dissociation constant for APC:axin binding	$K_{7}$	$50$	$n M$
Dissociation constant for β-catenin:DC binding	$K_{8}$	$120$	$n M$
Phosphorylation rate of β-catenin	$k_{9}$	$206$	$m i n^{- 1}$
Rate of phosphorylated β-catenin release from DC	$k_{10}$	$206$	$m i n^{- 1}$
Degradation rate of phosphorylated β-catenin	$k_{11}$	$0.417$	$m i n^{- 1}$
Synthesis rate of β-catenin	$v_{12}$	$0.423$	$n M m i n^{- 1}$
Degradation rate of β-catenin	$k_{13}$	$2.57 \cdot 10^{- 4}$	$m i n^{- 1}$
Synthesis rate of axin	$v_{14}$	$8.22 \cdot 10^{- 5}$	$n M m i n^{- 1}$
Degradation rate of axin	$k_{15}$	$0.167$	$m i n^{- 1}$
Dissociation constant for β-catenin:TCF binding	$K_{16}$	$30$	$n M$
Dissociation constant for β-catenin:APC binding	$K_{17}$	$1200$	$n M$
Total concentration of Disheveled	${D v l}_{t o t}$	$100$	$n M$
Total concentration of adenomatous polyposis coli	$A P C_{t o t}$	$100$	$n M$
Total concentration of T-cell factor	${T C F}_{t o t}$	$15$	$n M$
Total concentration of glycogen synthase kinase 3β	$G S K 3_{t o t}$	$50$	$n M$
Independent Variable		Label
Active Disheveled		$X_{2}$
APC/axin/GSK3 (* denotes phosphorylated)		$X_{3}$
APC/axin/GSK3		$X_{4}$
β-catenin/APC/axin*/GSK3		$X_{9}$
β-catenin*		$X_{10}$
β-catenin		$X_{11} (β c a t)$
axin		$X_{12}$
Dependent Variable		Label
Inactive Disheveled		$X_{1}$
GSK3		$X_{5}$
APC/axin		$X_{6}$
APC		$X_{7}$
β-catenin/APC/axin/GSK3		$X_{8}$
TCF		$X_{13}$
β-catenin/TCF		$X_{14}$
β-catenin/APC		$X_{15}$
Differential Equations
$\dot{[X_{2}]} = k_{1} \cdot W n t \cdot (D v l_{t o t} - [X_{2}]) - k_{2} \cdot [X_{2}]$
$(1 + \frac{[X_{11}]}{K_{8}}) \cdot \dot{[X_{3}]} + \frac{[X_{3}]}{K_{8}} \cdot \dot{[X_{11}]} = k_{4} \cdot [X_{4}] - k_{5} \cdot [X_{3}] - \frac{k_{9} \cdot [X_{3}] \cdot [X_{1} 1]}{K_{8}} + k_{10} \cdot [X_{9}]$
$\dot{[X_{4}]} = - (k_{3} \cdot [X_{2}] + k_{4} + k_{- 6}) \cdot [X_{4}] + k_{5} \cdot [X_{3}] + k_{6} \cdot G S K 3_{t o t} \cdot \frac{K_{17} \cdot [X_{12}] \cdot A P C_{t o t}}{K_{7} \cdot (K_{17} + [X_{11}])}$
$\dot{[X_{9}]} = \frac{k_{9} \cdot [X_{3}] \cdot [X_{11}]}{K_{8}} - k_{10} \cdot [X_{9}]$
$\dot{[X_{10}]} = k_{10} \cdot [X_{9}] - k_{11} \cdot [X_{10}]$
$(1 + \frac{[X_{3}]}{K_{8}} + \frac{K_{16} \cdot T C F_{t o t}}{(K_{16} + [X_{11}])^{2}} + \frac{K_{17} \cdot A P C_{t o t}}{(K_{17} + [X_{11}])^{2}}) \cdot \dot{[X_{11}]} + \frac{[X_{11}]}{K_{8}} \cdot \dot{[X_{3}]} = v_{12} - (\frac{k_{9} \cdot [X_{3}]}{K_{8}} + k_{13}) \cdot [X_{11}]$
$\begin{array}{ll} (1 + \frac{K_{17} \cdot A P C_{t o t}}{K_{7} \cdot (K_{17} + [X_{11}])}) \cdot \dot{[X_{12}]} - \frac{K_{17} \cdot [X_{12}] \cdot A P C_{t o t}}{K_{7} \cdot (K_{17} + [X_{11}])^{2}} \cdot \dot{[X_{11}]} \\ = k_{3} \cdot [X_{2}] \cdot [X_{4}] - k_{6} \cdot G S K_{t o t} \cdot \frac{K_{17} \cdot [X_{12} \cdot A P C_{t o t}]}{K_{7} \cdot (K_{17} + [X_{11}])} + k_{- 6} \cdot [X_{4}] + v_{14} \\ - k_{15} \cdot [X_{12}] \end{array}$
Equations for fast equilibrium reactions
$[X_{1}] = D v l_{t o t} - [X_{2}]$
$[X_{5}] = G S K 3_{t o t}$
$[X_{6}] = \frac{K_{17} \cdot [X_{12} \cdot A P C_{t o t}]}{K_{7} \cdot (K_{17} + [X_{11}])}$
$[X_{7}] = \frac{K_{17} \cdot A P C_{t o t}}{K_{17} \cdot (K_{17} + [X_{11}])}$
$[X_{8}] = \frac{[X_{3}] \cdot [X_{11}]}{K_{8}}$
$[X_{13}] = \frac{K_{16} \cdot T C F_{t o t}}{K_{16} + [X_{11}]}$
$[X_{14}] = \frac{[X_{11}] \cdot T C F_{t o t}}{K_{16} + [X_{11}]}$
$[X_{15}] = \frac{[X_{11}] \cdot A P C_{t o t}}{K_{17} + [X_{11}]}$

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

Parameter	Label	Value
Forward rate for Raf:RasGTP binding	$k_{3}$	$1.67 \cdot 10^{- 6}$	$m o l e c u l e^{- 1} s^{- 1}$
Reverse rate for Raf:RasGTP binding	$k_{b 3}$	$5.3 \cdot 10^{- 3}$	$s^{- 1}$
Phosphorylation rate for Raf by RasGTP	$k_{4}$	$1$	$s^{- 1}$
Forward rate of pRaf:P1 binding	$k_{7}$	$1.18 \cdot 10^{- 4}$	${m o l e c u l e}^{- 1} s^{- 1}$
Reverse rate of pRaf:P1 binding	$k_{b 7}$	$0.2$	$s^{- 1}$
Dephosphorylation rate of pRaf by P1	$k_{8}$	$1$	$s^{- 1}$
Forward rate of MEK:pRaf binding	$k_{9}$	$1.95 \cdot 10^{- 5}$	${m o l e c u l e}^{- 1} s^{- 1}$
Reverse rate of MEK:pRaf binding	$k_{b 9}$	$3.3 \cdot 10^{- 2}$	$s^{- 1}$
Phosphorylation rate of MEK by pRaf	$k_{10}$	$3.5$	$s^{- 1}$
Forward rate of pMEK:pRaf binding	$k_{11}$	$1.95 \cdot 10^{- 5}$	${m o l e c u l e}^{- 1} s^{- 1}$
Reverse rate of pMEK:pRaf binding	$k_{b 11}$	$3.3 \cdot 10^{- 2}$	$s^{- 1}$
Phosphorylation rate of pMEK by pRaf	$k_{12}$	$2.9$	$s^{- 1}$
Forward rate of dpMEK:P2 binding	$k_{13}$	$2.38 \cdot 10^{- 5}$	${m o l e c u l e}^{- 1} s^{- 1}$
Reverse rate of dpMEK:P2 binding	$k_{b 13}$	$0.8$	$s^{- 1}$
Dephosphorylation rate of dpMEK by P2	$k_{14}$	$5.8 \cdot 10^{- 2}$	$s^{- 1}$
Forward rate of pMEK:P2 binding	$k_{15}$	$4.5 \cdot 10^{- 7}$	${m o l e c u l e}^{- 1} s^{- 1}$
Reverse rate of pMEK:P2 binding	$k_{b 15}$	$0.5$	$s^{- 1}$
Dephosphorylation rate of pMEK by P2	$k_{16}$	$5.8 \cdot 10^{- 2}$	$s^{- 1}$
Forward rate of ERK:dpMEK binding	$k_{17}$	$8.9 \cdot 10^{- 5}$	${m o l e c u l e}^{- 1} s^{- 1}$
Reverse rate of ERK:dpMEK binding	$k_{b 17}$	$1.83 \cdot 10^{- 2}$	$s^{- 1}$
Phosphorylation rate of ERK by dpMEK	$k_{18}$	$16$	$s^{- 1}$
Forward rate of pERK:dpMEK binding	$k_{19}$	$8.9 \cdot 10^{- 5}$	${m o l e c u l e}^{- 1} s^{- 1}$
Reverse rate of pERK:dpMEK binding	$k_{b 19}$	$1.83 \cdot 10^{- 2}$	$s^{- 1}$
Phosphorylation rate of pERK by dpMEK	$k_{20}$	$5.7$	$s^{- 1}$
Forward rate of pERK:P3 binding	$k_{21}$	$8.33 \cdot 10^{- 6}$	${m o l e c u l e}^{- 1} s^{- 1}$
Reverse rate of pERK:P3 binding	$k_{b 21}$	$0.5$	$s^{- 1}$
Dephosphorylation rate of pERK by P3	$k_{22}$	$0.246$	$s^{- 1}$
Forward rate of dpERK:P3 binding	$k_{23}$	$2.35 \cdot 10^{- 5}$	${m o l e c u l e}^{- 1} s^{- 1}$
Reverse rate of dpERK:P3 binding	$k_{b 23}$	$0.6$	$s^{- 1}$
Dephosphorylation rate of dpERK by P3	$k_{24}$	$0.246$	$s^{- 1}$
Forward rate of Raf:dpERK binding	$k_{25}$	$1 \cdot 10^{- 6}$	${m o l e c u l e}^{- 1} s^{- 1}$
Reverse rate of Raf:dpERK binding	$k_{b 25}$	$1$	$s^{- 1}$
Hyper-phosphorylation rate of Raf by ppERK	$k_{26}$	$10$	$s^{- 1}$
Forward rate of pRaf:dpERK binding	$k_{27}$	$0$	${m o l e c u l e}^{- 1} s^{- 1}$
Reverse rate of pRaf:dpERK binding	$k_{b 27}$	$1$	$s^{- 1}$
Hyper-phosphorylation rate of phosphorylated Raf by dpERK	$k_{28}$	$10$	$s^{- 1}$
Forward rate of Rafi:P4 binding	$k_{29}$	$5 \cdot 10^{- 5}$	${m o l e c u l e}^{- 1} s^{- 1}$
Reverse rate of Rafi:P4 binding	$k_{b 29}$	$0.2$	$s^{- 1}$
Dephosphorylation rate of Rafi by P4	$k_{30}$	$0.5$	$s^{- 1}$
Total Raf	$R a f_{t o t}$	$4 \cdot 10^{4}$	$m o l e c u l e s$
Total MEK	$M E K_{t o t}$	$2.1 \cdot 10^{7}$	$m o l e c u l e s$
Total ERK	$E R K_{t o t}$	$2.21 \cdot 10^{7}$	$m o l e c u l e s$
Total phosphatase P1	$P 1_{t o t}$	$4 \cdot 10^{4}$	$m o l e c u l e s$
Total phosphatase P2	$P 2_{t o t}$	$4 \cdot 10^{5}$	$m o l e c u l e s$
Total phosphatase P3	$P 3_{t o t}$	$1 \cdot 10^{7}$	$m o l e c u l e s$
Total phosphatase P4	$P 4_{t o t}$	$4 \cdot 10^{4}$	$m o l e c u l e s$
Variable		Label
Unphosphorylated Raf		$R a f$
Raf bound to RasGTP		$R a f : R a s G T P$
Phosphorylated Raf		$p R a f$
Phosphatase for phosphorylated Raf		$P 1$
Phosphorylated Raf bound to its phosphatase		$p R a f : P 1$
Unphosphorylated MEK		$M E K$
MEK bound to its kinase		$M E K : p R a f$
Phosphorylated MEK		$p M E K$
Phosphorylated MEK bound to its kinase		$p M E K : p R a f$
Doubly-phosphorylated MEK		$d p M E K$
MEK phosphatase		$P 2$
Doubly-phosphorylated MEK bound to its phosphatase		$d p M E K : P 2$
Phosphorylated MEK bound to its phosphatase		$p M E K : P 2$
Unphosphorylated ERK		$E R K$
ERK bound to its kinase		$E R K : d p M E K$
Phosphorylated ERK		$p E R K$
Phosphorylated ERK bound to its kinase		$p E R K : d p M E K$
Doubly-phosphorylated ERK		$d p E R K$
ERK phosphatase		$P 3$
Phosphorylated ERK bound to its phosphatase		$p E R K : P 3$
Doubly-phosphorylated ERK bound to its phosphatase		$d p E R K : P 3$
Raf bound to doubly-phosphorylated ERK		$R a f : d p E R K$
Hyper-phosphorylated, ‘inactive’ Raf		$R a f i$
Phosphorylated Raf bound to doubly-phosphorylated ERK		$p R a f : d p E R K$
Phosphatase for hyper-phosphorylated Raf		$P 4$
Hyper-phosphorylated Raf bound to its phosphatase		$R a f i : P 4$
Differential Equations
$\begin{array}{ll} \dot{[R a f]} = - k_{3} \cdot [R a f] \cdot u (E G F) + k_{b 3} \cdot [R a f : R a s G T P] + k_{8} \cdot [p R a f : P 1] - k_{2} 5 \cdot [R a f] \\ \cdot [d p E R K] + k_{b 25} \cdot [R a f : d p E R K] + k_{30} \cdot [R a f i : P 4] \end{array}$
$\dot{[R a f : R a s G P T]} = k_{3} \cdot [R a f] \cdot u (E G F) - (k_{b 3} + k_{4}) \cdot [R a f : R a s G T P]$
$\begin{array}{ll} \dot{[p R a f]} = & k_{4} \cdot [R a f : R a s G T P] - k_{7} \cdot [p R a f] \cdot [P 1] + k_{b 7} \cdot [p R a f : P 1] - k_{9} \cdot [M E K] \cdot [p R a f] \\ + (k_{b 9} + k_{10}) \cdot [M E K : p R a f] - k_{11} \cdot [p M E K] \cdot [p R a f] + (k_{b 11} + k_{12}) \\ \cdot [p M E K : p R a f] - k_{27} \cdot [p R a f] \cdot [d p E R K] + k_{b 27} \cdot [p R a f : d p E R K] \end{array}$
$\dot{[P 1]} = - k_{7} \cdot [p R a f] \cdot [P 1] + (k_{b 7} + k_{8}) \cdot [p R a f : P 1]$
$\dot{[p R a f : P 1]} = k_{7} \cdot [p R a f] \cdot [P 1] - (k_{b 7} + k_{8}) \cdot [p R a f : P 1]$
$\dot{[M E K]} = - k_{9} \cdot [M E K] \cdot [p R a f] + k_{b 9} \cdot [M E K : p R a f] + k_{16} \cdot [p M E K : P 2]$
$\dot{[M E K : p R a f]} = k_{9} \cdot [M E K] \cdot [p R a f] - (k_{b 9} + k_{10}) \cdot [M E K : p R a f]$
$\begin{array}{ll} \dot{[p M E K]} = & k_{10} \cdot [M E K : p R a f] - k_{11} \cdot [p M E K] \cdot [p R a f] + k_{b 11} \cdot [p M E K : p R a f] + k_{14} \\ \cdot [d p M E K : P 2] - k_{15} \cdot [p M E K] \cdot [P 2] + k_{b 15} \cdot [p M E K : P 2] \end{array}$
$\dot{[p M E K : p R a f]} = k_{11} \cdot [p M E K] \cdot [p R a f] - (k_{b 11} + k_{12}) \cdot [p M E K : p R a f]$
$\begin{array}{ll} \dot{[d p M E K]} = & k_{12} \cdot [p M E K : p R a f] - k_{13} \cdot [d p M E K] \cdot [P 2] + k_{b 13} \cdot [d p M E K : P 2] - k_{17} \cdot [E R K] \\ \cdot [d p M E K] + (k_{b 17} + k_{18}) \cdot [E R K : d p M E K] - k_{19} \cdot [p E R K] \cdot [d p M E K] \\ + (k_{b 19} + k_{20}) \cdot [p E R K : d p M E K] \end{array}$
$\begin{array}{ll} \dot{[P 2]} = & - k_{13} \cdot [d p M E K] \cdot [P 2] + (k_{b 13} + k_{14}) \cdot [d p M E K : P 2] - k_{15} \cdot [p M E K] \cdot [P 2] \\ + (k_{b 15} + k_{16}) \cdot [p M E K : P 2] \end{array}$
$\dot{[d p M E K : P 2]} = k_{13} \cdot [d p M E K] \cdot [P 2] - (k_{b 13} + k_{14}) \cdot [d p M E K : P 2]$
$\dot{[p M E K : P 2]} = k_{15} \cdot [p M E K] \cdot [P 2] - (k_{b 15} + k_{16}) \cdot [p M E K : P 2]$
$\dot{[E R K]} = - k_{17} \cdot [E R K] \cdot [d p M E K] + k_{b 17} \cdot [E R K : d p M E K] + k_{22} \cdot [p E R K : P 3]$
$\dot{[E R K : d p M E K]} = k_{17} \cdot [E R K] \cdot [d p M E K] - (k_{b 17} + k_{18}) \cdot [E R K : d p M E K]$
$\begin{array}{ll} \dot{[p E R K]} = & k_{18} \cdot [E R K : d p M E K] - k_{19} \cdot [p E R K] \cdot [d p M E K] + k_{b 19} \cdot [p E R K : d p M E K] - k_{21} \\ \cdot [p E R K] \cdot [P 3] + k_{b 21} \cdot [p E R K : P 3] + k_{24} \cdot [d p E R K : P 3] \end{array}$
$\dot{[p E R K : d p M E K]} = k_{19} \cdot [p E R K] \cdot [d p M E K] - (k_{b 19} + k_{20}) \cdot [p E R K : d p M E K]$
$\begin{array}{ll} \dot{[d p E R K]} = & k_{20} \cdot [p E R K : d p M E K] - k_{23} \cdot [d p E R K] \cdot [P 3] + k_{b 23} \cdot [d p E R K : P 3] - k_{25} \cdot [R a f] \\ \cdot [d p E R K] + (k_{b 25} + k_{26}) \cdot [R a f : d p E R K] - k_{27} \cdot [p R a f] \cdot [d p E R K] \\ + (k_{b 27} + k_{28}) \cdot [p R a f : d p E R K] \end{array}$
$\begin{array}{ll} \dot{[P 3]} = & - k_{21} \cdot [p E R K] \cdot [P 3] + (k_{b 21} + k_{22}) \cdot [p E R K : P 3] - k_{23} \cdot [d p E R K] \cdot [P 3] \\ + (k_{b 23} + k_{24}) \cdot [d p E R K : P 3] \end{array}$
$\dot{[p E R K : P 3]} = k_{21} [p E R K] \cdot [P 3] - (k_{b 21} + k_{22}) \cdot [p E R K : P 3]$
$\dot{[d p E R K : P 3]} = k_{23} \cdot [d p E R K] \cdot [P 3] - (k_{b 23} + k_{24}) \cdot [d p E R K : P 3]$
$\dot{[R a f : d p E R K]} = k_{25} \cdot [R a f] \cdot [d p E R K] - (k_{b 25} + k_{26}) \cdot [R a f : d p E R K]$
$\dot{[R a f i]} = k_{26} \cdot [R a f : d p E R K] + k_{28} \cdot [p R a f : d p E R K] - k_{29} \cdot [R a f i] \cdot [P 4] + k_{b 29} \cdot [R a f i : P 4]$
$\dot{[p R a f : d p E R K]} = k_{27} \cdot [p R a f] \cdot [d p E R K] - (k_{b 27} + k_{28}) \cdot [p R a f : d p E R K]$
$\dot{[P 4]} = - k_{29} \cdot [R a f i] \cdot [P 4] + (k_{b 29} + k_{30}) \cdot [R a f i : P 4]$
$\dot{[R a f i : P 4]} = k_{29} \cdot [R a f i] \cdot [P 4] - (k_{b 29} + k_{30}) \cdot [R a f i : P 4]$
Algebraic Equations for conserved species
$R a f_{t o t} = [R a f] + [R a f : R a s G T P] + [p R a f] + [p R a f : P 1] + [M E K : p R a f] + [p M E K : p R a f] + [R a f : d p E R K] + [R a f i] + [p R a f : d p E R K] + [R a f i : P 4]$
$M E K_{t o t} = [M E K] + [M E K : p R a f] + [p M E K] + [p M E K : p R a f] + [d p M E K] + [d p M E K : P 2] + [p M E K : P 2] + [E R K : d p M E K] + [p E R K : d p M E K]$
$E R K_{t o t} = [E R K] + [E R K : d p M E K] + [p E R K] + [p E R K : d p M E K] + [d p E R K] + [p E R K : P 3] + [d p E R K : P 3] + [R a f : d p E R K] + [p R a f : d p E R K]$
$P 1_{t o t} = [P 1] + [p R a f : P 1]$
$P 2_{t o t} = [P 2] + [d p M E K : P 2] + [p M E K : P 2]$
$P 3_{t o t} = [P 3] + [p E R K : P 3] + [d p E R K : P 3]$
$P 4_{t o t} = [P 4] + [R a f i : P 4]$

Appendix 1—table 3

Parameters, variables and equations of the Tgfβ model.

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

Parameter	Label	Value
Phosphorylation rate of Smad2	$k_{p h o s}$	$4.0 \cdot 10^{- 4}$	$n M^{- 1} s^{- 1}$
Dephosphorylation rate of Smad2	$k_{d e p h o s}$	$6.6 \cdot 10^{- 3}$	$n M^{- 1} s^{- 1}$
Nuclear import rate of Smad2	$k_{i n 2}$	$2.6 \cdot 10^{- 3}$	$s^{- 1}$
Nuclear export rate of Smad2	$k_{e x 2}$	$5.6 \cdot 10^{- 3}$	$s^{- 1}$
Nuclear import rate of Smad4	$k_{i n 4}$	$2.6 \cdot 10^{- 3}$	$s^{- 1}$
Nuclear export rate of Smad4	$k_{e x 4}$	$2.6 \cdot 10^{- 3}$	$s^{- 1}$
Smad complex import factor	$C I F$	$5.7$
Forward rate for Smad complex binding	$k_{o n}$	$1.8 \cdot 10^{- 3}$	$n M^{- 1} s^{- 1}$
Reverse rate for Smad complex binding	$k_{o f f}$	$1.6 \cdot 10^{- 2}$	$s^{- 1}$
Cytoplasmic to nuclear volume ratio	$a$	$2.3$
Total Smad2 (initialized to cytoplasm)	$S 2_{t o t}$	$73.0$	$n M$
Total Smad4 (initialized to cytoplasm)	$S 4_{t o t}$	$73.0$	$n M$
Total phosphatase in nucleus	$P P a s e$	$1$	$n M$
Total Receptors	$R_{t o t}$	$1$	$n M$
Variable		Label
Cytoplasmic Smad2		$S 2 c$
Cytoplasmic phosphorylated Smad2		$p S 2 c$
Cytoplasmic Smad4		$S 4 c$
Cytoplasmic Smad2:Smad4 complex		$S 24 c$
Cytoplasmic Smad2:Smad2 complex		$S 22 c$
Nuclear Smad2		$S 2 n$
Nuclear phosphorylated Smad2		$p S 2 n$
Nuclear Smad4		$S 4 n$
Nuclear Smad2:Smad4 complex		$S 24 n$
Nuclear Smad2:Smad2 complex		$S 22 n$
Differential Equations
$\dot{[S 2 c]} = - k_{p h o s} \cdot u (T g f β) \cdot [S 2 c] - k_{i n 2} \cdot [S 2 c] + k_{e x 2} \cdot [S 2 n]$
$\begin{array}{ll} \dot{[p S 2 c]} = & k_{p h o s} \cdot u (T f g β) \cdot [S 2 c] - k_{i n 2} \cdot [p S 2 c] - k_{o n} \cdot [p S 2 c] \cdot ([S 4 c] + 2 \cdot [p S 2 c]) + k_{o f f} \\ \cdot ([S 24 c] + 2 \cdot [S 22 c]) + k_{e x 2} \cdot [p S 2 n] \end{array}$
$\dot{[S 4 c]} = - k_{i n 4} \cdot [S 4 c] - k_{o n} \cdot [p S 2 c] \cdot [S 4 c] + k_{o f f} \cdot [S 24 c] + k_{e x 4} \cdot [S 4 n]$
$\dot{[S 24 c]} = k_{o n} [p S 2 c] \cdot [S 4 c] - k_{o f f} \cdot [S 24 c] - k_{i n 2} \cdot C I F \cdot [S 24 c]$
$\dot{[S 22 c]} = k_{o n} \cdot [p S 2 c] 2 - k_{o f f} \cdot [S 22 c] - k_{i n 2} \cdot C I F \cdot [S 22 c]$
$\dot{[S 2 n]} = a \cdot k_{i n 2} \cdot [S 2 c] - a \cdot k_{e x 2} \cdot [S 2 n] + k_{d e p h o s} \cdot P P a s e \cdot [p S 2 n]$
$\begin{array}{ll} \dot{[p S 2 n]} = & a \cdot k_{i n 2} \cdot [p S 2 c] - a \cdot k_{e x 2} \cdot [p S 2 n] - k_{d e p h o s} \cdot P P a s e \cdot [p S 2 n] - k_{o n} \cdot [p S 2 n] \\ \cdot ([S 4 n] + 2 \cdot [p S 2 n]) + k_{o f f} \cdot ([S 24 n] + 2 \cdot [S 22 n]) \end{array}$
$\dot{[S 4 n]} = a \cdot k_{i n 4} \cdot [S 4 c] - a \cdot k_{e x 4} \cdot [S 4 n] - k_{o n} \cdot [p S 2 n] \cdot [S 4 n] + k_{o f f} \cdot [S 24 n]$
$\dot{[S 24 n]} = a \cdot k_{i n 2} \cdot C I F \cdot [S 24 c] + k_{o n} \cdot [p S 2 n] \cdot [S 4 n] - k_{o f f} \cdot [S 24 n]$
$\dot{[S 22 n]} = a \cdot k_{i n 2} \cdot C I F \cdot [S 22 c] + k_{o n} \cdot [p S 2 n] 2 - k_{o f f} \cdot [S 22 n]$
Algebraic Equations for conserved species
$S 2_{t o t} = [S 2 c] + [p S 2 c] + [S 24 c] + 2 \cdot [S 22 c] + (2 \cdot [S 22 n] + [S 24 n] + [p S 2 n] + [S 2 n])$
$S 4_{t o t} = [S 4 c] + [S 24 c] + \frac{1}{a} ([S 24 n] + [S 4 n])$

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

Pathway	Systems where graded response has been observed	References
Live imaging of single cells
Tgfβ pathway	Mouse myoblasts	Frick et al., 2017; Warmflash et al. (2012)
	Human epidermal keratinocytes	Nicolás et al. (2004); Warmflash et al. (2012) ; Schmierer et al. (2008); Vizán et al. (2013)
	Human cervical epithelial cells	Nicolás et al. (2004)
	Human breast epithelial cells	Strasen et al. (2018)
Canonical Wnt pathway	Human embryonic kidney cells	Kafri et al. (2016) (this is the only published live single-cell imaging study in the Wnt pathway so far)
ERK pathway	Mouse fibroblasts	Toettcher et al. (2013)
	Mouse embryonic fibroblasts	Mackeigan et al. (2005)
	Human non-small cell lung carcinoma	Cheong et al., 2011
	Human mammary gland cells	Selimkhanov et al. (2014); Perrett et al., 2013 Perrett et al., 2013
	Human cervical epithelial cells	Voliotis et al. (2014); Whitehurst et al. (2004); Perrett et al., 2013 Perrett et al., 2013
	Human foreskin fibroblasts	Whitehurst et al. (2004)
Immunofluorescence and FACS studies
Tgfβ pathway	Xenopus embryo	Schohl and Fagotto (2002)
	Mouse testes	Itman et al. (2009)
	Zebrafish embryo	Dubrulle et al. (2015)
Canonical Wnt pathway	Xenopus embryo	Schneider et al. (1996); Fagotto and Gumbiner (1994); Schohl and Fagotto (2002)
	Mouse embyo	Aulehla et al., 2008
	Planaria	Sureda-Gómez et al., 2016
	Sea anemone embryo	Wikramanayake et al., 2003
ERK pathway	Chick embryo	Delfini et al. (2005)
	Xenopus embryo	Schohl and Fagotto (2002)
	Human T lymphocyte cells	Lin et al. (2009)
	Rat adrenal gland cells	Santos et al., 2007

Data availability

All data generated or analysed during this study are included in the manuscript and supporting files. Source data files have been provided for Figures 2 and 3.

References

1. Ahmed S
2. Grant KG
3. Edwards LE
4. Rahman A
5. Cirit M
6. Goshe MB
7. Haugh JM
(2014) Data-driven modeling reconciles kinetics of ERK phosphorylation, localization, and activity states
Molecular Systems Biology 10:718.
https://doi.org/10.1002/msb.134708
- PubMed
- Google Scholar
1. Andrews SS
2. Peria WJ
3. Yu RC
4. Colman-Lerner A
5. Brent R
(2016) Push-Pull and feedback mechanisms can align signaling system outputs with inputs
Cell Systems 3:444–455.
https://doi.org/10.1016/j.cels.2016.10.002
- PubMed
- Google Scholar
(2012) Dynamic regulation of Tgf-B signaling by Tif1γ: a computational approach
PLOS ONE 7:e33761.
https://doi.org/10.1371/journal.pone.0033761
- PubMed
- Google Scholar
1. Ashall L
2. Horton CA
3. Nelson DE
4. Paszek P
5. Harper CV
6. Sillitoe K
7. Ryan S
8. Spiller DG
9. Unitt JF
10. Broomhead DS
11. Kell DB
12. Rand DA
13. Sée V
14. White MR
(2009) Pulsatile stimulation determines timing and specificity of NF-κB-dependent transcription
Science 324:242–246.
https://doi.org/10.1126/science.1164860
- PubMed
- Google Scholar
1. Aulehla A
2. Wiegraebe W
3. Baubet V
4. Wahl MB
5. Deng C
6. Taketo M
7. Lewandoski M
8. Pourquié O
(2008) A beta-catenin gradient links the clock and wavefront systems in mouse embryo segmentation
Nature Cell Biology 10:186–193.
https://doi.org/10.1038/ncb1679
- PubMed
- Google Scholar
1. Avraham R
2. Yarden Y
(2011) Feedback regulation of EGFR signalling: decision making by early and delayed loops
Nature Reviews Molecular Cell Biology 12:104–117.
https://doi.org/10.1038/nrm3048
- PubMed
- Google Scholar
1. Becker V
2. Schilling M
3. Bachmann J
4. Baumann U
5. Raue A
6. Maiwald T
7. Timmer J
8. Klingmüller U
(2010) Covering a broad dynamic range: information processing at the erythropoietin receptor
Science 328:1404–1408.
https://doi.org/10.1126/science.1184913
- PubMed
- Google Scholar
1. Cheong R
2. Rhee A
3. Wang CJ
4. Nemenman I
5. Levchenko A
(2011) Information transduction capacity of noisy biochemical signaling networks
Science 334:354–358.
https://doi.org/10.1126/science.1204553
- PubMed
- Google Scholar
(2003) The deaf and the dumb: the biology of ErbB-2 and ErbB-3
Experimental Cell Research 284:54–65.
https://doi.org/10.1016/S0014-4827(02)00101-5
- PubMed
- Google Scholar
1. Clevers H
2. Nusse R
(2012) Wnt/β-catenin signaling and disease
Cell 149:1192–1205.
https://doi.org/10.1016/j.cell.2012.05.012
- PubMed
- Google Scholar
1. Cohen-Saidon C
2. Cohen AA
3. Sigal A
4. Liron Y
5. Alon U
(2009) Dynamics and variability of ERK2 response to EGF in individual living cells
Molecular Cell 36:885–893.
https://doi.org/10.1016/j.molcel.2009.11.025
- PubMed
- Google Scholar
(2016) Control theory meets synthetic biology
Journal of The Royal Society Interface 13:20160380.
https://doi.org/10.1098/rsif.2016.0380
- PubMed
- Google Scholar
(2005) Control of the segmentation process by graded MAPK/ERK activation in the chick embryo
PNAS 102:11343–11348.
https://doi.org/10.1073/pnas.0502933102
- PubMed
- Google Scholar
1. Dougherty MK
2. Müller J
3. Ritt DA
4. Zhou M
5. Zhou XZ
6. Copeland TD
7. Conrads TP
8. Veenstra TD
9. Lu KP
10. Morrison DK
(2005) Regulation of Raf-1 by direct feedback phosphorylation
Molecular Cell 17:215–224.
https://doi.org/10.1016/j.molcel.2004.11.055
- PubMed
- Google Scholar
(2015) Response to Nodal morphogen gradient is determined by the kinetics of target gene induction
eLife 4:e05042.
https://doi.org/10.7554/eLife.05042
- PubMed
- Google Scholar
1. Fagotto F
2. Gumbiner BM
(1994)
Beta-catenin localization during Xenopus embryogenesis: accumulation at tissue and somite boundaries

Development 120:3667–3679.
- PubMed
- Google Scholar
1. Feliu E
2. Wiuf C
(2012) Variable elimination in chemical reaction networks with Mass-Action kinetics
SIAM Journal on Applied Mathematics 72:959–981.
https://doi.org/10.1137/110847305
- Google Scholar
1. Ferrell JE
2. Bhatt RR
(1997) Mechanistic studies of the dual phosphorylation of mitogen-activated protein kinase
Journal of Biological Chemistry 272:19008–19016.
https://doi.org/10.1074/jbc.272.30.19008
- PubMed
- Google Scholar
1. Ferrell JE
2. Machleder EM
(1998) The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes
Science 280:895–898.
https://doi.org/10.1126/science.280.5365.895
- PubMed
- Google Scholar
1. Frick CL
2. Yarka C
3. Nunns H
4. Goentoro L
(2017) Sensing relative signal in the Tgf-β/Smad pathway
PNAS 114:E2975–E2982.
https://doi.org/10.1073/pnas.1611428114
- PubMed
- Google Scholar
1. Fritsche-Guenther R
2. Witzel F
3. Sieber A
4. Herr R
5. Schmidt N
6. Braun S
7. Brummer T
8. Sers C
9. Blüthgen N
(2011) Strong negative feedback from Erk to Raf confers robustness to MAPK signalling
Molecular Systems Biology 7:489.
https://doi.org/10.1038/msb.2011.27
- PubMed
- Google Scholar
1. Gao S
2. Alarcón C
3. Sapkota G
4. Rahman S
5. Chen PY
6. Goerner N
7. Macias MJ
8. Erdjument-Bromage H
9. Tempst P
10. Massagué J
(2009) Ubiquitin ligase Nedd4L targets activated Smad2/3 to limit TGF-β signaling
Molecular Cell 36:457–468.
https://doi.org/10.1016/j.molcel.2009.09.043
- PubMed
- Google Scholar
1. Gerhart J
(1999) 1998 warkany lecture: signaling pathways in development
Teratology 60:226–239.
https://doi.org/10.1002/(SICI)1096-9926(199910)60:4<226::AID-TERA7>3.0.CO;2-W
- PubMed
- Google Scholar
1. Goentoro L
2. Kirschner MW
(2009) Evidence that fold-change, and not absolute level, of β-catenin dictates Wnt signaling
Molecular Cell 36:872–884.
https://doi.org/10.1016/j.molcel.2009.11.017
- PubMed
- Google Scholar
(2011) Studying Smad2 intranuclear diffusion dynamics by mathematical modelling of FRAP experiments
Integrative Biology 3:197–207.
https://doi.org/10.1039/c0ib00098a
- PubMed
- Google Scholar
1. Hekman M
2. Fischer A
3. Wennogle LP
4. Wang YK
5. Campbell SL
6. Rapp UR
(2005) Novel C-Raf phosphorylation sites: serine 296 and 301 participate in Raf regulation
FEBS Letters 579:464–468.
https://doi.org/10.1016/j.febslet.2004.11.105
- PubMed
- Google Scholar
(2012) Kinetic responses of β-catenin specify the sites of Wnt control
Science 338:1337–1340.
https://doi.org/10.1126/science.1228734
- PubMed
- Google Scholar
(2002) The IκB-NF-κB signaling module: temporal control and selective gene activation
Science 298:1241–1245.
https://doi.org/10.1126/science.1071914
- PubMed
- Google Scholar
Book
1. Hoppler SP
2. Moon RT
(2014) Wnt Signaling in Development and Disease: Molecular Mechanisms and Biological Functions
John Wiley & Sons.
https://doi.org/10.1002/9781118444122
- Google Scholar
1. Huang CY
2. Ferrell JE
(1996) Ultrasensitivity in the mitogen-activated protein kinase cascade
PNAS 93:10078–10083.
https://doi.org/10.1073/pnas.93.19.10078
- PubMed
- Google Scholar
(2002) Nucleocytoplasmic shuttling of Smads 2, 3, and 4 permits sensing of TGF-beta receptor activity
Molecular Cell 10:283–294.
https://doi.org/10.1016/S1097-2765(02)00585-3
- PubMed
- Google Scholar
1. Itman C
2. Small C
3. Griswold M
4. Nagaraja AK
5. Matzuk MM
6. Brown CW
7. Jans DA
8. Loveland KL
(2009) Developmentally regulated SMAD2 and SMAD3 utilization directs activin signaling outcomes
Developmental Dynamics 238:1688–1700.
https://doi.org/10.1002/dvdy.21995
- PubMed
- Google Scholar
1. Kafri P
2. Hasenson SE
3. Kanter I
4. Sheinberger J
5. Kinor N
6. Yunger S
7. Shav-Tal Y
(2016) Quantifying β-catenin subcellular dynamics and cyclin D1 mRNA transcription during Wnt signaling in single living cells
eLife 5:e16748.
https://doi.org/10.7554/eLife.16748
- PubMed
- Google Scholar
1. Kimelman D
2. Xu W
(2006) β-Catenin destruction complex: insights and questions from a structural perspective
Oncogene 25:7482–7491.
https://doi.org/10.1038/sj.onc.1210055
- PubMed
- Google Scholar
(1984)
Relationship between epidermal growth factor receptor occupancy and mitogenic response. quantitative analysis using a steady state model system.

The Journal of Biological Chemistry 259:5623–5654.
- PubMed
- Google Scholar
1. Kolch W
(2005) Coordinating ERK/MAPK signalling through scaffolds and inhibitors
Nature Reviews Molecular Cell Biology 6:827–837.
https://doi.org/10.1038/nrm1743
- PubMed
- Google Scholar
(1996) Partnership between DPC4 and SMAD proteins in TGF-β signalling pathways
Nature 383:832–836.
https://doi.org/10.1038/383832a0
- PubMed
- Google Scholar
(2016) Negative feedback regulation of the ERK1/2 MAPK pathway
Cellular and Molecular Life Sciences 73:4397–4413.
https://doi.org/10.1007/s00018-016-2297-8
- PubMed
- Google Scholar
1. Lee E
2. Salic A
3. Krüger R
4. Heinrich R
5. Kirschner MW
(2003) The roles of APC and Axin derived from experimental and theoretical analysis of the Wnt pathway
PLOS Biology 1:e10.
https://doi.org/10.1371/journal.pbio.0000010
- PubMed
- Google Scholar
1. Lee RE
2. Walker SR
3. Savery K
4. Frank DA
5. Gaudet S
(2014) Fold change of nuclear NF-κB determines TNF-induced transcription in single cells
Molecular cell 53:867–879.
https://doi.org/10.1016/j.molcel.2014.01.026
- PubMed
- Google Scholar
1. Lin J
2. Harding A
3. Giurisato E
4. Shaw AS
(2009) KSR1 modulates the sensitivity of mitogen-activated protein kinase pathway activation in T cells without altering fundamental system outputs
Molecular and Cellular Biology 29:2082–2091.
https://doi.org/10.1128/MCB.01634-08
- PubMed
- Google Scholar
(2009) Wnt/β-catenin signaling: components, mechanisms, and diseases
Developmental Cell 17:9–26.
https://doi.org/10.1016/j.devcel.2009.06.016
- PubMed
- Google Scholar
(2005) Graded mitogen-activated protein kinase activity precedes switch-like c-Fos induction in mammalian cells
Molecular and Cellular Biology 25:4676–4682.
https://doi.org/10.1128/MCB.25.11.4676-4682.2005
- PubMed
- Google Scholar
(2005) Smad transcription factors
Genes & Development 19:2783–2810.
https://doi.org/10.1101/gad.1350705
- PubMed
- Google Scholar
1. Nevozhay D
2. Adams RM
3. Murphy KF
4. Josic K
5. Balázsi G
(2009) Negative autoregulation linearizes the dose-response and suppresses the heterogeneity of gene expression
PNAS 106:5123–5128.
https://doi.org/10.1073/pnas.0809901106
- PubMed
- Google Scholar
(2013) Transferring a synthetic gene circuit from yeast to mammalian cells
Nature Communications 4:1451.
https://doi.org/10.1038/ncomms2471
- PubMed
- Google Scholar
(2004) Analysis of Smad nucleocytoplasmic shuttling in living cells
Journal of Cell Science 117:4113–4125.
https://doi.org/10.1242/jcs.01289
- PubMed
- Google Scholar
1. Nusse R
2. Clevers H
(2017) Wnt/β-Catenin Signaling, Disease, and Emerging Therapeutic Modalities
Cell 169:985–999.
https://doi.org/10.1016/j.cell.2017.05.016
- PubMed
- Google Scholar
1. Olsman N
2. Goentoro L
(2016) Allosteric proteins as logarithmic sensors
Proceedings of the National Academy of Sciences 113:E4423–E4430.
https://doi.org/10.1073/pnas.1601791113
- PubMed
- Google Scholar
(2014) The EGFR demonstrates linear signal transmission
Integr. Biol. 6:736–742.
https://doi.org/10.1039/C4IB00062E
- PubMed
- Google Scholar
(2013) Signaling to ERK from ErbB1 and PKC: Feedback, Heterogeneity and Gating
The Journal of Biological Chemistry jbc. M113. 455345.
https://doi.org/10.1074/jbc.M113.455345
- Google Scholar
1. Pires-daSilva A
2. Sommer RJ
(2003) The evolution of signalling pathways in animal development
Nature Reviews Genetics 4:39–49.
https://doi.org/10.1038/nrg977
- PubMed
- Google Scholar
1. Rotherham M
2. El Haj AJ
(2015) Remote activation of the Wnt/β-catenin signalling pathway using functionalised magnetic particles
PLOS ONE 10:e0121761.
https://doi.org/10.1371/journal.pone.0121761
- PubMed
- Google Scholar
1. Saito-Diaz K
2. Chen TW
3. Wang X
4. Thorne CA
5. Wallace HA
6. Page-McCaw A
7. Lee E
(2013) The way Wnt works: components and mechanism
Growth Factors 31:1–31.
https://doi.org/10.3109/08977194.2012.752737
- PubMed
- Google Scholar
(2007) Growth factor-induced MAPK network topology shapes Erk response determining PC-12 cell fate
Nature Cell Biology 9:324–330.
https://doi.org/10.1038/ncb1543
- PubMed
- Google Scholar
1. Schmierer B
2. Hill CS
(2005) Kinetic analysis of Smad nucleocytoplasmic shuttling reveals a mechanism for transforming growth factor β-dependent nuclear accumulation of Smads
Molecular and Cellular Biology 25:9845–9858.
https://doi.org/10.1128/MCB.25.22.9845-9858.2005
- PubMed
- Google Scholar
(2008) Mathematical modeling identifies Smad nucleocytoplasmic shuttling as a dynamic signal-interpreting system
PNAS 105:6608–6613.
https://doi.org/10.1073/pnas.0710134105
- PubMed
- Google Scholar
(1996) β-catenin translocation into nuclei demarcates the dorsalizing centers in frog and fish embryos
Mechanisms of Development 57:191–198.
https://doi.org/10.1016/0925-4773(96)00546-1
- PubMed
- Google Scholar
(2002) Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors
Nature Biotechnology 20:370–375.
https://doi.org/10.1038/nbt0402-370
- PubMed
- Google Scholar
1. Schohl A
2. Fagotto F
(2002)
β-catenin, MAPK and Smad signaling during early Xenopus development

Development 129:37–52.
- PubMed
- Google Scholar
1. Selimkhanov J
2. Taylor B
3. Yao J
4. Pilko A
5. Albeck J
6. Hoffmann A
7. Tsimring L
8. Wollman R
(2014) Systems biology. Accurate information transmission through dynamic biochemical signaling networks
Science 346:1370–1373.
https://doi.org/10.1126/science.1254933
- PubMed
- Google Scholar
1. Spiess AN
2. Neumeyer N
(2010) An evaluation of R2 as an inadequate measure for nonlinear models in pharmacological and biochemical research: a Monte Carlo approach
BMC Pharmacology 10:6.
https://doi.org/10.1186/1471-2210-10-6
- PubMed
- Google Scholar
(2010) Cis-interactions between Notch and Delta generate mutually exclusive signalling states
Nature 465:86–90.
https://doi.org/10.1038/nature08959
- PubMed
- Google Scholar
1. Stamos JL
2. Weis WI
(2013) The β-catenin destruction complex
Cold Spring Harbor Perspectives in Biology 5:a007898.
https://doi.org/10.1101/cshperspect.a007898
- PubMed
- Google Scholar
1. Strasen J
2. Sarma U
3. Jentsch M
4. Bohn S
5. Sheng C
6. Horbelt D
7. Knaus P
8. Legewie S
9. Loewer A
(2018) Cell-specific responses to the cytokine TGFβ are determined by variability in protein levels
Molecular Systems Biology 14:e7733.
https://doi.org/10.15252/msb.20177733
- PubMed
- Google Scholar
1. Sturm OE
2. Orton R
3. Grindlay J
4. Birtwistle M
5. Vyshemirsky V
6. Gilbert D
7. Calder M
8. Pitt A
9. Kholodenko B
10. Kolch W
(2010) The mammalian MAPK/ERK pathway exhibits properties of a negative feedback amplifier
Science Signaling 3:ra90.
https://doi.org/10.1126/scisignal.2001212
- PubMed
- Google Scholar
(2016) Localization of planarian β-CATENIN-1 reveals multiple roles during anterior-posterior regeneration and organogenesis
Development 143:4149–4160.
https://doi.org/10.1242/dev.135152
- PubMed
- Google Scholar
1. Tamai K
2. Zeng X
3. Liu C
4. Zhang X
5. Harada Y
6. Chang Z
7. He X
(2004) A mechanism for Wnt coreceptor activation
Molecular Cell 13:149–156.
https://doi.org/10.1016/S1097-2765(03)00484-2
- PubMed
- Google Scholar
1. Thurley K
2. Tovey SC
3. Moenke G
4. Prince VL
5. Meena A
6. Thomas AP
7. Skupin A
8. Taylor CW
9. Falcke M
(2014) Reliable encoding of stimulus intensities within random sequences of intracellular Ca2+ spikes
Science Signaling 7:ra59.
https://doi.org/10.1126/scisignal.2005237
- PubMed
- Google Scholar
(2013) Using optogenetics to interrogate the dynamic control of signal transmission by the Ras/Erk module
Cell 155:1422–1434.
https://doi.org/10.1016/j.cell.2013.11.004
- PubMed
- Google Scholar
1. Vizán P
2. Miller DS
3. Gori I
4. Das D
5. Schmierer B
6. Hill CS
(2013) Controlling long-term signaling: receptor dynamics determine attenuation and refractory behavior of the TGF-β pathway
Science Signaling 6:ra106.
https://doi.org/10.1126/scisignal.2004416
- PubMed
- Google Scholar
(2014) Information transfer by leaky, heterogeneous, protein kinase signaling systems
PNAS 111:E326–E333.
https://doi.org/10.1073/pnas.1314446111
- PubMed
- Google Scholar
1. Wang J
2. Tucker-Kellogg L
3. Ng IC
4. Jia R
5. Thiagarajan PS
6. White JK
7. Yu H
(2014) The self-limiting dynamics of TGF-β signaling in silico and in vitro, with negative feedback through PPM1A upregulation
PLOS Computational Biology 10:e1003573.
https://doi.org/10.1371/journal.pcbi.1003573
- PubMed
- Google Scholar
1. Warmflash A
2. Zhang Q
3. Sorre B
4. Vonica A
5. Siggia ED
6. Brivanlou AH
(2012) Dynamics of TGF- signaling reveal adaptive and pulsatile behaviors reflected in the nuclear localization of transcription factor Smad4
PNAS 109:E1947–E1956.
https://doi.org/10.1073/pnas.1207607109
- Google Scholar
(2004) Stimulus-coupled spatial restriction of extracellular signal-regulated kinase 1/2 activity contributes to the specificity of signal-response pathways
Molecular and Cellular Biology 24:10145–10150.
https://doi.org/10.1128/MCB.24.23.10145-10150.2004
- PubMed
- Google Scholar
1. Wikramanayake AH
2. Hong M
3. Lee PN
4. Pang K
5. Byrum CA
6. Bince JM
7. Xu R
8. Martindale MQ
(2003) An ancient role for nuclear β-catenin in the evolution of axial polarity and germ layer segregation
Nature 426:446–450.
https://doi.org/10.1038/nature02113
- PubMed
- Google Scholar
1. Xu L
2. Massagué J
(2004) Nucleocytoplasmic shuttling of signal transducers
Nature Reviews Molecular Cell Biology 5:209–219.
https://doi.org/10.1038/nrm1331
- PubMed
- Google Scholar
1. Yoon S
2. Seger R
(2006) The extracellular signal-regulated kinase: multiple substrates regulate diverse cellular functions
Growth Factors 24:21–44.
https://doi.org/10.1080/02699050500284218
- PubMed
- Google Scholar
1. Yu RC
2. Pesce CG
3. Colman-Lerner A
4. Lok L
5. Pincus D
6. Serra E
7. Holl M
8. Benjamin K
9. Gordon A
10. Brent R
(2008) Negative feedback that improves information transmission in yeast signalling
Nature 456:755–761.
https://doi.org/10.1038/nature07513
- PubMed
- Google Scholar
1. Zeng X
2. Tamai K
3. Doble B
4. Li S
5. Huang H
6. Habas R
7. Okamura H
8. Woodgett J
9. He X
(2005) A dual-kinase mechanism for Wnt co-receptor phosphorylation and activation
Nature 438:873–877.
https://doi.org/10.1038/nature04185
- PubMed
- Google Scholar
1. Zeng X
2. Huang H
3. Tamai K
4. Zhang X
5. Harada Y
6. Yokota C
7. Almeida K
8. Wang J
9. Doble B
10. Woodgett J
11. Wynshaw-Boris A
12. Hsieh JC
13. He X
(2008) Initiation of Wnt signaling: control of Wnt coreceptor Lrp6 phosphorylation/activation via frizzled, dishevelled and axin functions
Development 135:367–375.
https://doi.org/10.1242/dev.013540
- PubMed
- Google Scholar

Decision letter

Wenying Shou

Reviewing Editor; Fred Hutchinson Cancer Research Center, United States
Aviv Regev

Senior Editor; Broad Institute of MIT and Harvard, United States
Wenying Shou

Reviewer; Fred Hutchinson Cancer Research Center, United States
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 (k₂₅) – 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

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

"This ORCID iD identifies the author of this article:" 0000-0002-9669-0039
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

"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

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

Reviewing Editor

Wenying Shou, Fred Hutchinson Cancer Research Center, United States

Reviewers

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

Version history

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

Copyright

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.