Signaling pathways as linear transmitters
 Cited 2
 Views 2,410
 Annotations
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 wellcharacterized molecular details allow us to relate the internal processes of each pathway to their inputoutput 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 inputoutput 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.001Introduction
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; PiresdaSilva 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, allornone 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 ligandreceptor 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 ligandreceptor 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 ligandreceptor activation, and the output as the level of transcriptional regulator.
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; SaitoDiaz 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. Ligandreceptor 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). Ligandreceptor 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 systemslevel 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 foldchange 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 theorybased variable elimination and dimensional analysis, to derive analytical or semianalytical solutions to the steadystate 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 posttranslational 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; FritscheGuenther et al., 2011]; Tgfβ refs. [Schmierer et al., 2008; GonzálezPé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 ligandreceptor 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 steadystate (Goentoro and Kirschner, 2009):
where the input function $\mathrm{u}=\mathrm{u}\left(\mathrm{W}\mathrm{n}\mathrm{t}\right)$ is the rate of inhibition of the destruction complex (DC) via Dishevelled/Dvl, a function of ligandreceptor 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. $\mathrm{\alpha}$ and $\mathrm{\gamma}$ in Equation 1 are dimensionless parameter groups defined in Equations 2 and 3: $\mathrm{\alpha}$ characterizes βcatenin degradation by the destruction complex, and $\mathrm{\gamma}$ 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 $\mathrm{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 $\mathrm{\alpha}~66$, $\mathrm{\gamma}~1.4$, and for maximal stimulation, $\mathrm{u}~6.0$. The large $\mathrm{\alpha}$ reflects how βcatenin stability is primarily dictated by the destruction complex, that is, $\mathrm{\alpha}/\mathrm{u}\gg 1$ means that nonAxindependent degradation is minimal, and $\mathrm{\alpha}/\mathrm{u}\gg \mathrm{\gamma}$ 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; SaitoDiaz et al., 2013; Hoppler and Moon, 2014). With $\mathrm{\alpha}/\mathrm{u}\gg 1+\mathrm{\gamma}$, Equation 1 simplifies to
with detailed derivations presented in Appendix 1. Therefore, within physiologically relevant parameter values, the steadystate βcatenin concentration becomes a linear function of the input $\mathrm{u}$ (red line, Figure 2D). The linear inputoutput relationship holds for the entire dynamic range of the model, until the system saturates at maximal stimulation ($\mathrm{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 $\mathrm{\alpha}$ is decreased, breaking the requirement $\mathrm{\alpha}/\mathrm{u}\gg 1+\mathrm{\gamma}$ (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), ligandreceptor input is transmitted via a cascade of protein phosphorylation (Kolch, 2005; Yoon and Seger, 2006). In particular, ligandreceptor 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). Doublyphosphorylated ERK (dpERK) is a transcriptional regulator that affects a broad array of genes (Yoon and Seger, 2006). The multistep 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; CohenSaidon 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 hyperphosphorylation 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; FritscheGuenther 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 steadystate in terms of core variables (described in Appendix 1). We derived an analytical relationship between the steadystate output of the pathway $\left[\mathrm{d}\mathrm{p}\mathrm{E}\mathrm{R}\mathrm{K}\right]}_{\mathrm{s}\mathrm{s}$ and the input to the phosphorylation cascade $\mathrm{u}$:
Detailed derivations of Equation 5 are presented in Appendix 1. The input $\mathrm{u}=\mathrm{u}\left(\mathrm{E}\mathrm{G}\mathrm{F}\right)$ in Equation 5 is the concentration of active Ras, which is activated via GTP loading at the ligandreceptor complex (Kolch, 2005). The parameter groups $\mathrm{\alpha}$, $\mathrm{\beta}$, $\mathrm{\gamma}$, and $\mathrm{\delta}$ in Equation 5 are defined in Equations 6–9, where the ellipses indicate additional small terms (expanded in Appendix 1). The relative magnitudes of $\mathrm{\alpha}$, $\mathrm{\beta}$, $\mathrm{\gamma}$, and $\mathrm{\delta}$ indicate how the Raf pool partitions during signaling (Equations A21, A29–A31). The dimensionless group $\mathrm{\alpha}\cdot \mathrm{u}$ relates to the amount of free, phosphorylated Raf ($\mathrm{\alpha}$, blueshaded in Figure 2B), $\mathrm{\beta}\cdot {\left[\mathrm{d}\mathrm{p}\mathrm{E}\mathrm{R}\mathrm{K}\right]}_{\mathrm{s}\mathrm{s}}$ describes the amount of Raf inhibited through negative feedback by dpERK ($\mathrm{\beta}$, redshaded in Figure 2B), $\mathrm{\delta}$ relates to the amount of unphosphorylated ($\mathrm{\delta}$, blueshaded in Figure 2B), and $\mathrm{\gamma}\cdot \mathrm{u}$ relates to the amount of phosphorylated Raf bound to other proteins (e.g. to MEK, brownshaded in Figure 2B). Equation 5 is not a closed solution, as it includes the term ${\left[\mathrm{p}\mathrm{R}\mathrm{a}\mathrm{f}\right]}_{\mathrm{s}\mathrm{s}}$, and there are variables included in parameter groups $\mathrm{\alpha}$, $\mathrm{\beta}$, $\mathrm{\gamma}$. 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 (FritscheGuenther et al., 2011; Dougherty et al., 2005). This means that $\mathrm{\beta}\cdot {\left[\mathrm{d}\mathrm{p}\mathrm{E}\mathrm{R}\mathrm{K}\right]}_{\mathrm{s}\mathrm{s}}~\left(\mathrm{\alpha}+\mathrm{\gamma}\right)\cdot \mathrm{u}+\mathrm{\delta}$. 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 $\mathrm{K}$ as the relative change of ${\left[\mathrm{d}\mathrm{p}\mathrm{E}\mathrm{R}\mathrm{K}\right]}_{\mathrm{s}\mathrm{s}}$ with respect to ${\left[\mathrm{p}\mathrm{R}\mathrm{a}\mathrm{f}\right]}_{\mathrm{s}\mathrm{s}}$, ultrasensitivity entails that $\mathrm{K}\gg 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 $\mathrm{\beta}{\cdot \left[\mathrm{d}\mathrm{p}\mathrm{E}\mathrm{R}\mathrm{K}\right]}_{\mathrm{s}\mathrm{s}}~\left(\mathrm{\alpha}+\mathrm{\gamma}\right)\cdot \mathrm{u}+\mathrm{\delta}$ and $\mathrm{K}\gg 1$, the negative feedback holds the level of pRaf constant (${\left[\mathrm{p}\mathrm{R}\mathrm{a}\mathrm{f}\right]}_{\mathrm{s}\mathrm{s}}\approx {\mathrm{R}}_{\mathrm{s}}$, details in Appendix 1). With these two features, strong negative feedback and ultrasensitivity, dpERK becomes a linear function of the input $\mathrm{u}$:
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 $\mathrm{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 ligandreceptor interactions is transmitted by the Smad proteins. There are several classes of Smad proteins, including the receptorregulated Smads (RSmads) and the common Smad (coSmad or Smad4) (Massagué et al., 2005). Ligandactivated receptors phosphorylate RSmads. Phosphorylated RSmads bind to the coSmad, and shuttle into the nucleus and regulate broad target genes. In the nucleus, the Smad complex dissociates and RSmads are constitutively dephosphorylated 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 RSmad2 data, the general architecture of signal transmission is conserved across all five RSmads (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 steadystate concentration of Smad complex in the nucleus:
In Equation 11, the input function $\mathrm{u}=\mathrm{u}\left(\mathrm{T}\mathrm{g}\mathrm{f}\mathrm{\beta}\right)$ is the active fraction of Tgfβ receptors. The parameter $\mathrm{a}$ is the nucleocytoplasmic volume ratio. The dimensionless parameter groups $\mathrm{\alpha}$, $\mathrm{\beta}$, and $\mathrm{\gamma}$ in Equation 11 are defined in Equations 12–14, where the ellipses indicate additional small terms (expanded in Appendix 1). $\mathrm{\alpha}$, $\mathrm{\beta}$, and $\mathrm{\gamma}$ describe how the Smad2 pool partitions during signaling (Equations A44, A50, A51): $\mathrm{\alpha}\cdot \mathrm{u}$ relates to the amount of nuclear Smad complex ($\mathrm{\alpha}$, blueshaded in Figure 2C, captures the parameters related to complex formation and translocation to the nucleus), $\mathrm{\beta}$ relates to the amount of free, unphosphorylated Smad2 ($\mathrm{\beta}$, redshaded in Figure 2C, captures the parameters related to complex dissociation and translocation to the cytoplasm), and $\mathrm{\gamma}\cdot \mathrm{u}$ loosely relates to the remaining Smad2 pool ($\mathrm{\gamma}$ is brownshaded in Figure 2C). Phosphorylated Smad2 quickly forms complex (Lagna et al., 1996), so $\mathrm{\beta}$ essentially corresponds to total monomeric Smad2. Finally, Equation 11 is not a closed solution, since variable ${\left[\mathrm{S}4\mathrm{n}\right]}_{\mathrm{s}\mathrm{s}}$ appears in $\mathrm{\alpha}$. We numerically tested that it is constant within 2% for nonsaturating 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 nonsaturating inputs ($\mathrm{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 $\mathrm{\beta}~46$, $\mathrm{\alpha}\cdot \mathrm{u}~1.5$, and $\mathrm{\gamma}\cdot \mathrm{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 ($\mathrm{\beta}\gg \left(\mathrm{\alpha}+\mathrm{\gamma}\right)\cdot \mathrm{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 $\mathrm{\beta}\gg \left(\mathrm{\alpha}+\mathrm{\gamma}\right)\cdot \mathrm{u}$, the first term in the denominator of Equation 11 is small, and concentration of nuclear Smad complex becomes a linear function of input:
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 RSmad phosphatase is inhibited such that $\mathrm{\beta}~\left(\mathrm{\alpha}+\mathrm{\gamma}\right)\cdot \mathrm{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 inputoutput 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 ligandreceptor 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 steadystate within 6 hr (Figure 3—figure supplement 1). Accordingly, all subsequent measurements were done at 6 hr after Wnt stimulation.
To measure the inputoutput 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).
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 MichaelisMenten kinetics that describe ligand binding in the model, we confirmed that the linearity does not extend upstream to Wnt dose: both phosphoLRP5/6 and βcatenin show nonlinear response to Wnt dose (Figure 3—figure supplement 2). Therefore, in the Wnt pathway, a nonlinear ligandreceptor processing step is followed by linear signal transmission through the core intracellular pathway.
Next, to measure the inputoutput 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 ligandreceptor processing. Here, we specifically probe linearity in the core transmission step of the pathway. Detecting the input level, EGFactivated 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 doublyphosphorylated 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 inputoutput 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 12fold 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 singlecell level (Figure 3—figure supplement 4), in agreement with a previous singlecell, live imaging study (CohenSaidon 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. $\mathrm{\alpha}/\mathrm{u}\gg 1+\mathrm{\gamma}$). This regime is not infinite: for instance, a tenfold decrease in $\mathrm{\alpha}$ (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 inputoutput relationship, βcatenin vs. phosphoLRP5/6 level, up to 90% of maximal phosphoLRP5/6 input (blue circles, Figure 3C). As expected, we found that inhibiting the destruction complex (decreasing $\mathrm{\alpha}$ in the model) reduced the range of linearity. The nontreated cells (blue circles, Figure 3A) exhibit a linear inputoutput relationship over a 4.4fold range of LRP input, whereas the CHIRtreated cells show a linear inputoutput relationship over only a 2.8fold 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 CHIRtreated 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 dualfunction 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 dualrole 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 inputoutput behavior in the Wnt pathway, as predicted by our analytics, and second, that GSK3β coregulation 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 $\mathrm{\beta}$ 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 Raf1 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 Raf1 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 Raf1 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 inputoutput 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 wellestablished 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).
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 doseresponse 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; CohenSaidon et al., 2009; Lee et al., 2014; Thurley et al., 2014; Frick et al., 2017). In linear inputoutput systems, the stimulated output correlates linearly to the basal output; thus, the foldchange 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 foldchange in the level of transcriptional regulator (Goentoro and Kirschner, 2009; CohenSaidon et al., 2009; Lee et al., 2014; Frick et al., 2017). Therefore, selecting for linearity may naturally confer the benefits of superposition, doseresponse alignment, and a robust foldchange 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 foldchange detection at the receptorlevel (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, allornone 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 protocolpBABEpuroCRAF that contains the wt human Raf1 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 sitedirected mutagenesis kit (New England Biolabs, E0554S). The mutant and wt Raf1 were then placed downstream of a CMV promoter.
Cell lines and cell culture
Request a detailed protocolRKO cells (ATCC, CRL2577) and H1299 cells (ATCC, CRL5803) 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 mMLglutamine (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 mMLglutamine (Invitrogen). Both cell lines tested negative for mycoplasma contamination.
Transfection of Raf1 constructs
Request a detailed protocolH1299 cells were transfected with the mutant and wt Raf1 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 protocolThe following antibodies were purchased from Cell Signaling Technologies: antiPhosphop44/42 MAPK (Erk1/2) (Thr202/Tyr204) (E10) Mouse mAb #9106, antihistone H3 (D1H2) XP Rabbit mAb #4499, anticRaf Antibody #9422, antiphosphoLRP6 (Ser1490) Antibody #2568, antiGAPDH (D4C6R) Mouse mAb #97166. AntiBetacatenin mouse mAb was purchased from BD Transduction Laboratories (#610153) and antiGAPDH rabbit antibody was purchased from Abcam (ab9485). The following fluorescent secondary antibodies were purchased from Fisher Scientific: IRDye 800CW Goat antiMouse IgG (926–32210) and IRDye 680LT Goat antiRabbit IgG (92668021).
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 protocolRKO cells were pretreated 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 protocolRKO cells at 70% confluency were scraped in PBS, pelleted, and snapfrozen, and then thawed in NP40 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% BisTris Plus Gel (Thermofisher, NW04120BOX). H1299 cells at 70% confluence were scraped in NP40 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% TrisGlycine Mini Gel (ThermoFisher, XP04200BOX).
Quantitative Western blots
Request a detailed protocolProteins 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 backgroundsubtracted 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 phosphoLRP5/6, variation in the foldactivation from experiment to experiment could artificially stretch the data along the x and yaxis, and introduce artifacts into the relationship between phosphoLRP5/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 inputoutput relationship: linearity was observed without normalization in cases where loading was already uniform (Figure 3—figure supplement 8).
L1 and L2norm analysis
Request a detailed protocolL1norm 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 L1norm is computed as the area between the polynomial fit and the diagonal. Linearity is defined in this context as L1norm < 0.1. L2norm analysis for Wnt pathway data was performed using a Pearson’s coefficient, and L2norm analysis for ERK pathway data was performed using the coefficient of correlation, $\mathrm{r}}^{2$.
Akaike information criterion
Request a detailed protocolTo score the validity of nonlinear model fits for Figure 3D, we used the biascorrected Akaike Information Criterion as described in ref. (Spiess and Neumeyer, 2010), which assesses goodnessoffit and model parsimony. The weighted Aikaike $\mathrm{w}\left(\mathrm{A}\mathrm{I}\mathrm{C}\right)$ 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 steadystates of the Tgfβ and ERK pathways. This technique was developed to handle the complexity of large chemical reaction networks. By eliminating variables from the steadystate solution, we can express the steadystate 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 steadystate solution consists of a large set of variables, each with a polynomial equation describing its steadystate.
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 firstorder 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.
In this network, there is a cut $\{\mathrm{A},\mathrm{B},\mathrm{C}\}$ that contains exactly one product and one reactant for each reaction. We have highlighted this cut in the reaction set:
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,
which are set to zero at steadystate, and two additional conservation equations:
The variable elimination technique allows us to reduce the steadystate 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 firstorder and homogenous with respect to our cut, and rewrite them in matrix form. We use the subscript ‘ss’ to denote steadystate:
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 firstorder 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 closedform analytical solutions for steadystate. The matrix formulation and variable elimination technique immediately provides us with a set of solvable variables. The solution to the matrix equation above is:
c is a scaling factor not constrained by the matrix equation. With the use of the conservation Equation 5, we can calculate $\mathrm{c}$ and express the steady state of all three species solely in terms of the parameters of the network, and ${\left[\mathrm{D}\right]}_{\mathrm{s}\mathrm{s}}$. For instance, the solution for $\left[\mathrm{C}\right]}_{\mathrm{s}\mathrm{s}$ is below.
The solutions for ${\left[\mathrm{A}\right]}_{\mathrm{s}\mathrm{s}}$, ${\left[\mathrm{B}\right]}_{\mathrm{s}\mathrm{s}}$, and ${\left[\mathrm{C}\right]}_{\mathrm{s}\mathrm{s}}$ derived from the variable elimination technique still depend on ${\left[\mathrm{D}\right]}_{\mathrm{s}\mathrm{s}}$. If we plug in the solutions for the cut species, we can obtain polynomial equations for the remaining species (in this case ${\left[\mathrm{D}\right]}_{\mathrm{s}\mathrm{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, ${\mathrm{k}}_{2}{\mathrm{k}}_{3}$, ${\mathrm{k}}_{1}{\mathrm{k}}_{3}{\left[\mathrm{D}\right]}_{\mathrm{s}\mathrm{s}}^{2}$, and ${\mathrm{k}}_{1}{\mathrm{k}}_{2}{\left[\mathrm{D}\right]}_{\mathrm{s}\mathrm{s}}^{2}$ represent the unnormalized fraction of ${\mathrm{T}}_{1}$ that exists as A, B, and C, respectively. The normalization factor for these fractions is $\mathrm{c}/{\mathrm{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 ${\mathrm{T}}_{1}$. For instance, increasing the value of ${\mathrm{k}}_{1}$ will increase the amount of ${\mathrm{T}}_{1}$ that exists as $\mathrm{B}$ and $\mathrm{C}$, while necessarily decreasing the amount of $\mathrm{A}$ (assuming ${\left[\mathrm{D}\right]}_{\mathrm{s}\mathrm{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 steadystate
We previously derived an expression for βcatenin in steadystate (Goentoro and Kirschner, 2009):
where the parameters are dimensionless groups of the binding rate constants and protein concentrations:
The input function $\mathrm{u}=\mathrm{u}\left(\mathrm{W}\mathrm{n}\mathrm{t}\right)$ corresponds to the rate at which Wnt stimulation inhibits the destruction complex, normalized by ${\mathrm{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:
Within the parameter regime measured in cells, the analytical expression for βcatenin dramatically simplifies. We can perform the following firstorder Taylor expansion:
This holds true for $\raisebox{1ex}{$\mathrm{\alpha}$}\!\left/ \!\raisebox{1ex}{$\mathrm{u}$}\right.\gg \mathrm{\gamma}$. Furthermore, we can make the approximation $1\mathrm{\gamma}+\raisebox{1ex}{$\mathrm{\alpha}$}\!\left/ \!\raisebox{1ex}{$\mathrm{u}$}\right.\approx \raisebox{1ex}{$\mathrm{\alpha}$}\!\left/ \!\raisebox{1ex}{$\mathrm{u}$}\right.$ as long as $\raisebox{1ex}{$\mathrm{\alpha}$}\!\left/ \!\raisebox{1ex}{$\mathrm{u}$}\right.\gg 1$ also holds. We can encompass these two inequalities within $\mathrm{\alpha}/\mathrm{u}\gg 1+\mathrm{\gamma}$. The equation simplifies to:
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. ${\mathrm{k}}_{25}$ characterizes the negative feedback from dpERK to unphosphorylated Raf, and ${\mathrm{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 ${\mathrm{k}}_{25}$, and set ${\mathrm{k}}_{27}$ to zero.
Solving the ERK model at steadystate
In the ERK pathway, doubly phosphorylated ERK is produced by the Raf/MEK/ERK cascade of phosphorylation,
There is a negative feedback within the pathway, such that,
where $\mathrm{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:
This allows us to express the steadystate 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:
With this, we derive the expression for ${\left[\mathrm{p}\mathrm{R}\mathrm{a}\mathrm{f}\right]}_{\mathrm{s}\mathrm{s}}$,
where the parameter groups are:
The ellipses indicate additional small terms (i.e. <10% of the previous terms, numerically calculated using the model parameters and $\mathrm{u}=4.5\mathrm{e}4$ molecules). All the calculations for this paper use these truncated parameter groups. The complete parameter groups are written below:
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, $\mathrm{\alpha}\cdot \mathrm{u}$ relates to the amount of free, phosphorylated Raf since $\mathrm{\alpha}\cdot \mathrm{u}/\left(\left(\mathrm{\alpha}+\mathrm{\gamma}\right)\cdot \mathrm{u}+\mathrm{\beta}{\left[\mathrm{d}\mathrm{p}\mathrm{E}\mathrm{R}\mathrm{K}\right]}_{\mathrm{s}\mathrm{s}}+\mathrm{\delta}\right)$ is the fraction of Raf present as pRaf. Thus, as $\mathrm{\alpha}\cdot \mathrm{u}$ increases relative to $\mathrm{\gamma}\cdot \mathrm{u}+\mathrm{\beta}{\left[\mathrm{d}\mathrm{p}\mathrm{E}\mathrm{R}\mathrm{K}\right]}_{\mathrm{s}\mathrm{s}}+\mathrm{\delta}$, the amount of pRaf also increases.
We can define three subpopulations of Raf: Raf inhibited by dpERK, $\left[{\mathrm{R}}_{\mathrm{i}}\right]$; Raf activated by RasGTP (input),$\left[{\mathrm{R}}_{\mathrm{a}}\right]$; and unphosphorylated Raf $\left[{\mathrm{R}}_{\mathrm{n}}\right]$. Specifically:
We can calculate the steadystate of each subpopulation as:
Thus, in the same sense that $\mathrm{\alpha}\cdot \mathrm{u}$ relates to the amount of free phosphorylated Raf, $\mathrm{\beta}\cdot {\left[\mathrm{d}\mathrm{p}\mathrm{E}\mathrm{R}\mathrm{K}\right]}_{\mathrm{s}\mathrm{s}}$relates to the amount of inhibited Raf, $\mathrm{\gamma}\cdot \mathrm{u}$ relates to the amount of phosphorylated Raf bound to other proteins (not free), and $\mathrm{\delta}$ 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 ${\left[\mathrm{d}\mathrm{p}\mathrm{E}\mathrm{R}\mathrm{K}\right]}_{\mathrm{s}\mathrm{s}}=\mathrm{g}\left({\left[\mathrm{p}\mathrm{R}\mathrm{a}\mathrm{f}\right]}_{\mathrm{s}\mathrm{s}}\right)$ 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.3fold change in pRaf leads to a 9fold change in $\mathrm{d}\mathrm{p}\mathrm{E}\mathrm{R}\mathrm{K}$ (from 10% to 90% of max, Figure 2—figure supplement 1BC).
We therefore approximate ${\left[\mathrm{p}\mathrm{R}\mathrm{a}\mathrm{f}\right]}_{\mathrm{s}\mathrm{s}}$ by a value ${\mathrm{R}}_{\mathrm{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 ${\left[\mathrm{d}\mathrm{p}\mathrm{E}\mathrm{R}\mathrm{K}\right]}_{\mathrm{s}\mathrm{s}}$ becomes a linear function of input:
Lastly, we write the value of two terms in Equation A32 below, numerically calculated using the parameter values of the model:
We can neglect the second term, yielding:
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
We can now derive a general expression for the relative change of ${\left[\mathrm{p}\mathrm{R}\mathrm{a}\mathrm{f}\right]}_{\mathrm{s}\mathrm{s}}$ with respect to a relative change in $\mathrm{u}$. We use the notation $\mathrm{d}\widehat{\mathrm{x}}=\mathrm{d}\mathrm{ln}\mathrm{x}=\mathrm{d}\mathrm{x}/\mathrm{x}$.
Next, we define the response coefficient $\mathrm{K}$ between ${\left[\mathrm{d}\mathrm{p}\mathrm{E}\mathrm{R}\mathrm{K}\right]}_{\mathrm{s}\mathrm{s}}$ and${\left[\mathrm{p}\mathrm{R}\mathrm{a}\mathrm{f}\right]}_{\mathrm{s}\mathrm{s}}$:
From Equation A21, we get the partial derivatives:
Using these two equations, we find that:
When $\mathrm{K}\gg 1$ and $\mathrm{\beta}{\left[\mathrm{d}\mathrm{p}\mathrm{E}\mathrm{R}\mathrm{K}\right]}_{\mathrm{s}\mathrm{s}}\mathrm{}~\mathrm{}\left(\mathrm{\alpha}+\mathrm{\gamma}\right)\cdot \mathrm{u}+\mathrm{\delta}$, we see that
Therefore, ${\left[\mathrm{p}\mathrm{R}\mathrm{a}\mathrm{f}\right]}_{\mathrm{s}\mathrm{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 ${\left[\mathrm{d}\mathrm{p}\mathrm{E}\mathrm{R}\mathrm{K}\right]}_{\mathrm{s}\mathrm{s}}$ becomes a linear function of input.
It is not guaranteed that the system is stable as $\mathrm{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 inputoutput linearity. In this model, induction of the output species $\mathrm{E}$ is a twostep process:
An input $\mathrm{u}$ increases the amount of species $\mathrm{R}$, which in turn influences $\mathrm{E}$ as $\mathrm{E}=\mathrm{g}\left(\mathrm{R}\right)$. There is negative feedback from $\mathrm{E}$ to $\mathrm{R}$, which in the limit of strong negative feedback is inversely proportional to $\mathrm{E}$.
Next, we specify the function $\mathrm{g}\left(\mathrm{R}\right)$ such that $\mathrm{K}={\mathrm{K}}_{0}$, where $\mathrm{K}$ is the relative change of $\mathrm{E}$ with respect to $\mathrm{R}$. As ${\mathrm{K}}_{0}$ increases, therefore, the function $\mathrm{g}\left(\mathrm{R}\right)$ becomes more ultrasensitive.
Solving for $\mathrm{E}$, we see that in the limit of $\mathrm{K}={\mathrm{K}}_{0}\gg 1$, $\mathrm{E}$ becomes a linear function of $\mathrm{u}$, and $\mathrm{R}$ is held constant at R_{s}.
While we do not have an explicit function for $\mathrm{g}\left(\mathrm{R}\right)$ 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 $\mathrm{g}\left(\mathrm{R}\right)$ in the region where $\mathrm{K}\gg 1$. We also show that these results hold outside the limit of strong negative feedback, as long as the feedbackinhibited pool of R is comparable to the remaining pool.
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 steadystate
We use the variable elimination technique described in section ‘Variable Elimination’ to derive an analytical expression for the steadystate concentration of nuclear Smad complex. First, based on the measured parameter values, and as confirmed by simulations, the extent of Smad2Smad2 binding is limited. We therefore neglect this reaction in subsequent analysis. We identify the following cut of the Tgfβ model:
which is subject to the conservation equation:
Thus, we can eliminate these variables from the steadystate polynomial solution, with dependence only on variables outside this cut:
Using this relationship, we derive an expression for the nuclear Smad complex$(\mathrm{S}24\mathrm{n}$) at steadystate,
where the parameter groups are:
Here the input function $\mathrm{u}=\mathrm{u}\left(\mathrm{T}\mathrm{g}\mathrm{f}\mathrm{\beta}\right)$ 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 ${\left[\mathrm{S}4\mathrm{c}\right]}_{\mathrm{s}\mathrm{s}}$ and ${\left[\mathrm{S}4\mathrm{n}\right]}_{\mathrm{s}\mathrm{s}}$ calculated for $\mathrm{u}=0$). All calculations for the paper use these truncated parameter groups. The complete parameter groups are written below:
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, $\mathrm{\alpha}\cdot \mathrm{u}$ relates to the amount of nuclear Smad complex since $\mathrm{\alpha}\cdot \mathrm{u}/\left(\left(\mathrm{\alpha}+\mathrm{\gamma}\right)\cdot \mathrm{u}+\mathrm{\beta}\right)$ is the fraction of Smad2 present as S24n. Thus, as $\mathrm{\alpha}\cdot \mathrm{u}$ increases relative to$\mathrm{\gamma}\cdot \mathrm{u}+\mathrm{\beta}$, the amount of S24n also increases.
By definition, the parameter groups $\mathrm{\beta}$ and $\mathrm{\gamma}\cdot \mathrm{u}$ capture the remaining inputindependent and inputdependent polynomials, respectively. Nevertheless, we would like to understand the physical significance of the parameter groups. We can calculate the amount of unphosphorylated Smad2 as:
$\mathrm{\delta}$ captures the dependence of nuclear, unphosphorylated Smad on the input. With the measured parameters, $\mathrm{\beta}\gg \mathrm{\delta}\cdot \mathrm{u}$, so we have
This means that $\mathrm{\beta}$ relates to the amount of unphosphorylated Smad2 in the same sense that $\mathrm{\alpha}\cdot \mathrm{u}$ relates to nuclear Smad complex. We can also express the remaining Smad2 species as:
However, as $\mathrm{\delta}$ is of the same order of magnitude as $\mathrm{\gamma}$, the parameter group $\mathrm{\gamma}$ 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
where we have used a nonsaturating input ($\mathrm{u}=\mathrm{}0.2$). Therefore, with measured parameters, $\mathrm{\beta}\gg \left(\mathrm{\alpha}+\mathrm{\gamma}\right)\cdot \mathrm{u}$. With this, the denominator in the $[\mathrm{S}4{\mathrm{n}]}_{\mathrm{s}\mathrm{s}}$ equation simplifies, and the concentration of Smad complex becomes a linear function of the input:
References
 1
 2
 3
 4
 5

6
Feedback regulation of EGFR signalling: decision making by early and delayed loopsNature Reviews Molecular Cell Biology 12:104–117.https://doi.org/10.1038/nrm3048
 7
 8

9
The deaf and the dumb: the biology of ErbB2 and ErbB3Experimental Cell Research 284:54–65.https://doi.org/10.1016/S00144827(02)001015
 10
 11

12
Control theory meets synthetic biologyJournal of The Royal Society Interface 13:20160380.https://doi.org/10.1098/rsif.2016.0380
 13

14
Regulation of Raf1 by direct feedback phosphorylationMolecular Cell 17:215–224.https://doi.org/10.1016/j.molcel.2004.11.055
 15

16
Betacatenin localization during Xenopus embryogenesis: accumulation at tissue and somite boundariesDevelopment 120:3667–3679.

17
Variable elimination in chemical reaction networks with MassAction kineticsSIAM Journal on Applied Mathematics 72:959–981.https://doi.org/10.1137/110847305

18
Mechanistic studies of the dual phosphorylation of mitogenactivated protein kinaseJournal of Biological Chemistry 272:19008–19016.https://doi.org/10.1074/jbc.272.30.19008
 19
 20

21
Strong negative feedback from Erk to Raf confers robustness to MAPK signallingMolecular Systems Biology 7:489.https://doi.org/10.1038/msb.2011.27
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31

32
Developmentally regulated SMAD2 and SMAD3 utilization directs activin signaling outcomesDevelopmental Dynamics 238:1688–1700.https://doi.org/10.1002/dvdy.21995
 33
 34

35
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.

36
Coordinating ERK/MAPK signalling through scaffolds and inhibitorsNature Reviews Molecular Cell Biology 6:827–837.https://doi.org/10.1038/nrm1743
 37

38
Negative feedback regulation of the ERK1/2 MAPK pathwayCellular and Molecular Life Sciences 73:4397–4413.https://doi.org/10.1007/s0001801622978
 39
 40
 41
 42

43
Graded mitogenactivated protein kinase activity precedes switchlike cFos induction in mammalian cellsMolecular and Cellular Biology 25:4676–4682.https://doi.org/10.1128/MCB.25.11.46764682.2005
 44
 45

46
Transferring a synthetic gene circuit from yeast to mammalian cellsNature Communications 4:1451.https://doi.org/10.1038/ncomms2471

47
Analysis of Smad nucleocytoplasmic shuttling in living cellsJournal of Cell Science 117:4113–4125.https://doi.org/10.1242/jcs.01289
 48

49
Allosteric proteins as logarithmic sensorsProceedings of the National Academy of Sciences 113:E4423–E4430.https://doi.org/10.1073/pnas.1601791113

50
The EGFR demonstrates linear signal transmissionIntegr. Biol. 6:736–742.https://doi.org/10.1039/C4IB00062E

51
Signaling to ERK from ErbB1 and PKC: Feedback, Heterogeneity and GatingThe Journal of Biological Chemistry jbc. M113. 455345.https://doi.org/10.1074/jbc.M113.455345

52
The evolution of signalling pathways in animal developmentNature Reviews Genetics 4:39–49.https://doi.org/10.1038/nrg977
 53

54
The way Wnt works: components and mechanismGrowth Factors 31:1–31.https://doi.org/10.3109/08977194.2012.752737
 55
 56
 57

58
βcatenin translocation into nuclei demarcates the dorsalizing centers in frog and fish embryosMechanisms of Development 57:191–198.https://doi.org/10.1016/09254773(96)005461
 59

60
βcatenin, MAPK and Smad signaling during early Xenopus developmentDevelopment 129:37–52.
 61
 62
 63

64
The βcatenin destruction complexCold Spring Harbor Perspectives in Biology 5:a007898.https://doi.org/10.1101/cshperspect.a007898

65
Cellspecific responses to the cytokine TGFβ are determined by variability in protein levelsMolecular Systems Biology 14:e7733.https://doi.org/10.15252/msb.20177733
 66
 67

68
A mechanism for Wnt coreceptor activationMolecular Cell 13:149–156.https://doi.org/10.1016/S10972765(03)004842
 69
 70
 71
 72
 73
 74
 75
 76

77
Nucleocytoplasmic shuttling of signal transducersNature Reviews Molecular Cell Biology 5:209–219.https://doi.org/10.1038/nrm1331
 78
 79
 80
 81
Decision letter

Wenying ShouReviewing Editor; Fred Hutchinson Cancer Research Center, United States

Aviv RegevSenior Editor; Broad Institute of MIT and Harvard, United States

Wenying ShouReviewer; Fred Hutchinson Cancer Research Center, United States

Steven S AndrewsReviewer; 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 inputoutput 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 switchlike responses, how would those parameters and inputs fare in your analytics? Would they have told you that the response will be nonlinear?
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 inputoutput relations in several model signaling pathways, which were the canonical Wnt, ERK, and Tgfβ pathways. The modeling started with established ODEbased 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 xaxis 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 doseresponse 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 inputoutput 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 doseresponse 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 NFkappaB pathway is a good example. Why was the NFkappaB 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 inputoutput 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 L1norm, as described in PMID:19279212 and PMID:23385595, which should be cited. It would be necessary to apply the L1norm 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 rearrange 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 L1norm) 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 inputoutput 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 foldchange 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 L1norm). 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 wildtype (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.034Author 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 inputoutput 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 switchlike responses, how would those parameters and inputs fare in your analytics? Would they have told you that the response will be nonlinear?
Switchlike 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., L1 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 inputoutput relations in several model signaling pathways, which were the canonical Wnt, ERK, and Tgfβ pathways. The modeling started with established ODEbased 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 xaxis 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 Wntdependent 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 EGFactivated RasGTP. To clarify this, we have now made u(ligand) more obvious in two ways:
A) We include u(ligand) function in the xaxis of Figure 2DF 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 2AC 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 doseresponse 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 foldchange 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 systemslevel behaviors, they do not necessarily reproduce finetuned 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 finetuned aspects such as the response gain (i.e., slope of the inputoutput relationship).
Consistently, our aim in Figure 3 was to assess whether the measured inputoutput 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 adhoc nonlinear functions, we used the full ERK model. We first fitted the gain of the model to the data (i.e., the yrange), and afterward, the only parameter varied was the strength of the negative feedback (k_{25}) – 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 inputoutput 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 doseresponse 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 NFkappaB pathway is a good example. Why was the NFkappaB 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 inputoutput 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, 5051). Numerical simulations of a wellestablished 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 inputoutput 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 proteincoupled 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, BaradaranMashinchi and Torres, 2018), are not conserved in mammals. Finally, with regard to the specific feedback we are considering here, where ERK hyperphosphorylates 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 RafMEKERK 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 L1norm, as described in PMID:19279212 and PMID:23385595, which should be cited. It would be necessary to apply the L1norm 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 (L2norm). As there is complementarity between L1 and L2norm, we also follow reviewer 3’s suggestion and added L1norm 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 foldchange 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 inputoutput 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 L1norm). 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, foldchange 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 rearrange 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 L1norm) 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 inputoutput 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 L1norm < 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 inputoutput 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 singlecell dynamics is now increasingly done, double quantitation of inputoutput 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 singlecell 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 singlecell ERK response (as discussed in their paper); however, the dynamic range of the sensor was not large enough to distinguish between doseresponse models (<2fold, compared to 12fold for our Western Blot data). Another method, quantitative dualimmunofluorescence, 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, nonsmall cell lung carcinoma, have been observed to show graded response at singlecell 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 L1norm). 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 theoryexperiment 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 finetuned 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 (phophoLRP), regardless of uncertainty in externally measured Wnt dose, it is notable that phosphoLRP5/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 wildtype (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 04 ng/mL EGF). We now include in the main text a data point for 3 ng/mL EGF in wt cells as well. Performing L1norm analysis over 04 ng/mL EGF, we find for wt cells, L1norm = 0.03, and for mutant cells, L1norm = 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 (phosphoLRP 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.035Article and author information
Author details
Funding
James S. McDonnell Foundation (220020365)
 Lea Goentoro
National Science Foundation (NSF.145863)
 Lea Goentoro
National Institutes of Health (5T32GM00761637)
 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
Publication 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
© 2018, Nunns et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,410
 Page views

 570
 Downloads

 2
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Download citations (links to download the citations from this article in formats compatible with various reference manager tools)
Open citations (links to open the citations from this article in various online reference manager services)
Further reading

 Cell Biology
 Immunology and Inflammation

 Cell Biology