Signaling pathways as linear transmitters

Abstract

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

Keywords: cell biology; computational biology; human cell lines; linearity; signal transmission; signaling pathways; systems biology.

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

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

Figure 2.. The Wnt, ERK, and Tgfβ pathways are linear signal transmitters.

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

Figure 2—figure supplement 1.. Model simulations for the ERK pathway.

(A) Parameter groups in the ERK model are constant to within 10%, over the physiologically relevant range of u considered here, justifying the inclusion of variables into the parameter groups. (B–C) The dpERK output is an ultrasensitive function of both free and total phosphorylated Raf. The values ${\mathrm{E}}_{\mathrm{s}}$ and ${\mathrm{R}}_{\mathrm{s}}$ are illustrated in (B), and are defined in Appendix 1. (D–F) Numerical simulation of pulsatile response in the ERK pathway. (D) A pulse of input, RasGTP, is generated by EGF addition in an ERK model that includes details of receptor desensitization (Schoeberl et al., 2002). Basal activity of Ras is included to ensure constitutive negative feedback (Fritsche-Guenther et al., 2011). (E) dpERK output also exhibits a pulsatile response, peaking within 10 min. (F) We plot the peak dpERK output against peak input for a range of physiologically relevant u(EGF) doses, and find that it matches our steady-state predictions for linear input-output behavior. (G) Five-fold Raf overexpression does not break the linear input-output behavior.

Figure 2—figure supplement 2.. The predicted linearity extends throughout the dynamic range of the ERK and Tgfβ pathways.

(A–B) Numerical simulation of the ERK and Tgfβ models. (A) The ERK model shows linear input-output relationship up to 93% of dpERK activation. (B) The Tgfβ pathway shows linear input-output relationship throughout the entire input range (from 0 to 1). Linearity was analyzed using the L1-norm or least absolute deviation (see Materials and methods). The blue range indicates where L1-norm was computed.

Figure 2—figure supplement 3.. Model simulations for the Tgfβ pathway.

(A) Nuclear Smad4 concentration is constant to within 2%, over a physiologically relevant range of u(Tgfβ) considered here, justifying its inclusion into parameter group $\mathrm{\alpha }$. (B–D) Numerical simulation of pulsatile response in the Tgfβ pathway. (B) A pulse of input, active Tgfβ receptor, is generated by Tgfβ addition in a model that includes details of receptor desensitization (Vizán et al., 2013). (C) S24n output also exhibits a pulsatile response. (D) We plot the peak S24n output against peak input, and find that it matches our steady-state predictions for linear input-output behavior.

Figure 2—figure supplement 4.. Incorporating into the Wnt model the dual function of GSK3β in phosphorylating β-catenin and LRP5/6.

We include the role of GSK3β in phosphorylating LRP5/6 into the input function u(Wnt), such that u(Wnt) is a function of GSK3β and the phosphorylation rate k₉ (for simplicity, the same rate as GSK3β phosphorylation of β-catenin), and the reverse rate k_r. In both models, β-catenin increases in response to GSK3β inhibition (e.g., by CHIR99021). However, only the model with the dual function of GSK3β shows a decrease in input range that we observed experimentally.

Figure 2—figure supplement 5.. The requirements for linear signal transmission in the Wnt, Tgfβ, and ERK pathway.

In each plot, we varied S, defined in the equation shown on the x-axis, and simulated the input-output curve over the dynamic range of the model. The parameters in the equations are as defined in the main text. For the ERK and Tgfβ pathway, $\mathrm{\alpha }$ and γ are linked in such a way that they could not easily be varied independently. Linearity was assessed using the L-1 norm, which ranges from 0 to 0.5, with L-1 norm < 0.1 indicating linearity. L1-norm analysis was performed over the full dynamic range of the system, i.e., u(Wnt) = 1- 6, u(ERK) = 0 to 110,000 molecules of Ras-GTP, which gave 90% activation of [dpERK] in unperturbed cells, and u(Tgfβ) = 0 to 1.

Figure 2—figure supplement 6.. Linear signal transmission occurs over a range of parameters in the model.

In this analysis, the parameter groups in each model were varied as indicated e.g., 3x is three-fold increase, 0.3x is three-fold decrease. 1x corresponds to the measured parameters. Plotted in each box is the input-output relationship, numerically simulated over the full dynamic range of the models, i.e., 1-6 for u(Wnt), 0-10⁵ for u(EGF), and 0-1 for u(Tgfβ). For simplicity, all outputs are normalized from 0 to 1. Grey shade: the unperturbed state. Purple shade: linear input-output response, as defined by L-1 norm < 0.1.

Figure 2—figure supplement 7.. Numerical simulation of the input-output relationship of the NF-κB pathway.

We used the model first built by Hoffman et al. in 2002 (Hoffmann et al., 2002) and later revised by Ashall et al. in 2009 (Ashall et al., 2009). The parameters in the model have been measured or fitted to single-cell dynamics in multiple cell types (Hoffmann et al., 2002; Ashall et al., 2009). We simulated the model here over a physiologically observed dynamic range, i.e., Lee et al. (2014) observed in HeLa cells that at saturating ligand dose (10 ng/mL TNFα, set to one in the model),~25% of NF-κB pool is nuclear. Linearity is assessed using the L1-norm, where L1-norm < 0.1 indicates linear relationship (see Materials and methods).

Figure 3.. Linearity was observed experimentally in the Wnt and ERK pathways.

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

Figure 3—figure supplement 1.. LRP5/6 phosphorylation and β-catenin accumulation are already at steady state at 6 hr after Wnt stimulation.

RKO cells were treated with 160 ng/mL Wnt3A for the specified times, and then assayed for phospho-LRP5/6 and β-catenin level by Western blot. Error bars are standard error of the mean from 2 to 4 biological replicates. Data are plotted relative to the sample at time zero, and normalized to the average maximal activation across experiments. The grey lines connect the mean of each time point.

Figure 3—figure supplement 2.. The dynamic range of Wnt signaling in RKO cells.

RKO cells were treated with the specified dose of Wnt3A for six hours, and then assayed for phospho-LRP5/6 and β-catenin by quantitative Western blot. Data are plotted relative to unstimulated samples. (A–B) In wt cells, phospho-LRP5/6 (A) shows > 90% of maximal response at 200 ng/mL Wnt3A, while β-catenin (B) shows 70% of maximal response at 200 ng/mL, and subsequently incremental response until 640 ng/mL Wnt3A. (C–D) In cells pre-treated with 1 µM CHIR99021, phospho-LRP5/6 (C) shows > 90% of maximal response at 80 ng/mL Wnt3A, while β-catenin (D) shows 70% of maximal response at 80 ng/mL and continues incremental activation at higher doses. (E). Cells were treated with 160 ng/mL Wnt3A and assayed for phospho-LRP5/6. Cells pre-treated with 1 µM CHIR99021 (N = 3) exhibited 50% the level of phospho-LRP5/6 as untreated cells (N = 3). The grey lines simply connect the means of data.

Figure 3—figure supplement 3.. ERK activation peaks at 5 min after EGF stimulation.

H1299 cells were treated with 1 ng/mL EGF for the specified times, and then assayed for dpERK1/2 by Western blot. Data is plotted relative to the samples at time zero, with at least three biological replicates per time point.

Figure 3—figure supplement 4.. Single-cell immunofluorescence measurements show graded ERK response to EGF.

In these experiments, H1299 cells were treated with varying doses of EGF for 5 min, and then fixed and analyzed for immunofluorescence against doubly phosphorylated ERK (dpERK). (A) Representative images of cells treated with the indicated doses of EGF. (B) The intensity of nuclear level of dpERK staining across individual cells. Cell nuclei were delineated using DAPI staining (for EGF doses 0, 0.3, 1.3, and 2.0 ng/mL, N = 453, 381, 373, and 413 cells, respectively).

Figure 3—figure supplement 5.. WT Raf-1 overexpression does not affect linear dose-response.

H1299 cells over-expressing Raf-1 were treated with the indicated dose of EGF for 5 min, and then assayed for dpERK1/2 by Western blot. The grey line is a fit from a linear model with r² = 0.99. Data is plotted relative to unstimulated samples, with total N = 9.

Figure 3—figure supplement 6.. Expression of Raf S29/289/296/301/642A induces non-linear dose-response.

H1299 cells expressing the Raf mutant Raf S29/289/296/301/642 were treated with the indicated dose of EGF for 5 min, and then assayed for dpERK1/2 by Western blot. Data is plotted relative to unstimulated samples, with total N = 5.

Figure 3—figure supplement 7.. Technical variability from Western blot.

(A) The level of β-catenin and phosphorylated LRP, measured across different lanes. (B) Ligand-stimulated change in β-catenin and phophorylated LRP level, measured in six independent Western blots. CV is coefficient of variation, defined as standard deviation/mean.

Figure 3—figure supplement 8.. Linearity is not an artifact of loading control normalization.

In these two independent experiments, RKO cells were stimulated with a range of Wnt3A dose (0-160 ng/mL), the cells lysed after 6 hr, and analyzed for Western blot against β-catenin and phosphorylated LRP5/6 (pLRP5/6). Top row: In each experiment, GAPDH intensity varies with <10% CV across samples. Bottom row: Raw β-catenin and LRP intensity data without normalization with GAPDH loading control. The measurements are plotted relative to unstimulated cells. Grey lines are least squares regression lines, and $\mathrm{\rho }$ is the Pearson correlation coefficient.

Figure 3—figure supplement 9.. Linearity was observed across independent experiments.

(A) In these two independent experiments, RKO cells were stimulated with a range of Wnt3A doses, lysed after 6 hr, and analyzed for Western blot against β-catenin and phospho-LRP5/6. (B) In these four independent experiments, H1299 cells were stimulated with a range of EGF doses, lysed after 5 min, and analyzed for Western blot against doubly-phosphorylated ERK. All measurements are plotted relative to unstimulated cells. Grey lines are least squares regression, ρ is Pearson correlation coefficient, and r² is correlation coefficient.

Figure 4.. Benefits of linearity.

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

Appendix 1—scheme 1.. Network of four proteins.

Appendix 1—scheme 2.. Reaction set corresponding to protein network.

Appendix 1—scheme 3.. Toy model of the ERK pathway.