Research ArticleCell Biology

Context-specific flow through the MEK/ERK module produces cell- and ligand-specific patterns of ERK single and double phosphorylation

See allHide authors and affiliations

Sci. Signal.  02 Feb 2016:
Vol. 9, Issue 413, pp. ra13
DOI: 10.1126/scisignal.aab1967

ERK phosphorylation patterns

In the ERK pathway, the dual-specificity kinase MEK phosphorylates a threonine and a tyrosine residue in ERK, and this dual-phosphorylated form is the fully active kinase. Iwamoto et al. used mass spectrometry, quantitative Western blotting, and mathematical modeling to explore MEK-dependent phosphorylation dynamics of ERK in skin and liver cells exposed to either a cytokine, IL-6, or a growth factor, HGF. Not surprisingly, the different stimuli produced different dynamics of ERK phosphorylation, and skin and liver cells responded differently to the same ligand. The dynamics of the changes in the abundance of the phosphorylated forms of ERK (pT-ERK, pY-ERK, and pTpY-ERK) and the relative distributions of the single- and double-phosphorylated forms of ERK were different. Mathematical modeling indicated that distinct network structures with or without regulated feedback loops produced the different dynamics of ERK phosphorylation and distributions of phosphorylated ERK. This study provides biochemical insight into how a single pathway can produce distinct responses, such as differentiation or proliferation.


The same pathway, such as the mitogen-activated protein kinase (MAPK) pathway, can produce different cellular responses, depending on stimulus or cell type. We examined the phosphorylation dynamics of the MAPK kinase MEK and its targets extracellular signal–regulated kinase 1 and 2 (ERK1/2) in primary hepatocytes and the transformed keratinocyte cell line HaCaT A5 exposed to either hepatocyte growth factor or interleukin-6. By combining quantitative mass spectrometry with dynamic modeling, we elucidated network structures for the reversible threonine and tyrosine phosphorylation of ERK in both cell types. In addition to differences in the phosphorylation and dephosphorylation reactions, the HaCaT network model required two feedback mechanisms, which, as the experimental data suggested, involved the induction of the dual-specificity phosphatase DUSP6 and the scaffold paxillin. We assayed and modeled the accumulation of the double-phosphorylated and active form of ERK1/2, as well as the dynamics of the changes in the monophosphorylated forms of ERK1/2. Modeling the differences in the dynamics of the changes in the distributions of the phosphorylated forms of ERK1/2 suggested that different amounts of MEK activity triggered context-specific responses, with primary hepatocytes favoring the formation of double-phosphorylated ERK1/2 and HaCaT A5 cells that produce both the threonine-phosphorylated and the double-phosphorylated form. These differences in phosphorylation distributions explained the threshold, sensitivity, and saturation of the ERK response. We extended the findings of differential ERK phosphorylation profiles to five additional cultured cell systems and matched liver tumor and normal tissue, which revealed context-specific patterns of the various forms of phosphorylated ERK.


Cells respond to multiple external cues using receptors and a complex, highly interlinked intracellular signaling network. Understanding cellular information processing requires the comparative analysis of differences in the signaling networks that specifically produce distinct biological responses.

A biochemical modification that mediates signal propagation is the change in the phosphorylation status of signaling molecules. An example of this type of modification is the highly conserved mitogen-activated protein kinase (MAPK) pathway in which varying the timing, duration, and amplitude of activation controls fundamental cellular processes, such as proliferation and differentiation (1). Extracellular signal–regulated kinase 1 and 2 (collectively, ERK1/2) are key components of the three-tiered MAPK pathway, which relays extracellular signals from cell surface receptors, such as cytokine receptors and receptor tyrosine kinases.

An essential feature of this signaling cascade is protein activation by multisite phosphorylation. ERK is only active if doubly phosphorylated (25), a status that is accomplished by dual phosphorylation by the upstream MAPK kinase (MEK). Multisite phosphorylation is a widespread feature of signal transduction, coordinating the specificity and sensitivity of cell signaling (6). Multisite phosphorylation can result in a cycle if the phosphorylation and dephosphorylation reactions result in different phosphorylated intermediates. Multisite phosphorylation mechanisms can produce thresholds for activation and switch-like behavior, also referred to as ultrasensitivity (79). Dual phosphorylation of ERK by MEK may proceed in a processive or a distributive manner. A processive mechanism, which involves a single interaction between the substrate and the kinase that results in phosphorylation of all targeted residues, would achieve full activation of ERK within a very short time span and facilitate rapid signal transduction (10), whereas distributive phosphorylation, which requires multiple interactions between the substrate and kinase, would provide higher sensitivity to small changes in the concentration of active MEK. A distributive mechanism would also enable signal amplification and a kinetic proofreading mechanism (11).

Previous systems biology studies have examined the entire MAPK network (1220), and the effect of scaffolds (21, 22) and negative feedback (23, 24) have been identified. Scaffold proteins, such as paxillin (25), localize kinases and their substrates in close proximity, which potentially accelerates signal transmission and enables processive phosphorylation. Negative feedback regulators, including dual-specificity phosphatases (DUSPs) (26), contribute to the robustness of a signaling response. Here, we focused on the context-specific phosphorylation of ERK by MEK (the MEK/ERK module) in two experimental model systems, primary murine hepatocytes (PMHs) and the human benign keratinocyte cell line HaCaT-Ras A5 (HaCaT A5), upon exposure to either hepatocyte growth factor (HGF) or interleukin-6 (IL-6). Although cells from different tissues, these cells are both from tissues that can regenerate. However, an important difference is that PMHs are nontransformed, whereas HaCaT A5 cells are transformed by Ras, which results in a low basal activation of the ERK pathway. Both of these cells respond to HGF: PMHs proliferate (27, 28) and HaCaT A5 cells proliferate and migrate (29). Both cell types also respond to IL-6. In PMHs, IL-6 induces DNA synthesis (27) and has been reported to be essential for liver regeneration (30). In HaCaT A5 cells, IL-6 stimulates proliferation and migration (31). Thus, these model systems represent cells from two regenerative tissues (32, 33), one that is transformed with a constitutively active ERK pathway due to mutation and one that has a more typical basal state for the ERK pathway.

We performed comparative modeling analysis, which provided evidence for a distributive, cyclic ERK phosphorylation mechanism and the presence of two HaCaT-specific regulatory feedback events that were not present in PMHs. We found that the stimulation induced changes in the distribution of the ERK forms: The dynamic changes in the amounts of nonphosphorylated (np-ERK), monophosphorylated (pT-ERK or pY-ERK), and double-phosphorylated (pTpY-ERK) ERK differed between the two types of cells and in response to the two ligands. Whereas PMHs preferentially produced double-phosphorylated ERK in response to HGF or IL-6, the HaCaT A5 cells produced similar amounts of pT-ERK and pTpY-ERK. Dose-response analyses of the dependency of the distribution of ERK phosphorylated forms on active MEK revealed that cell type–specific network structures could explain the threshold, sensitivity, and saturation behavior of the MEK/ERK module in both cells. Finally, we extended the findings of differential ERK phosphorylation profiles to primary keratinocytes, additional keratinocyte cell lines, primary human hepatocytes, a human liver carcinoma cell line, and matched liver tumor and normal tissue to gain insight into the context-specific patterns of the various forms of phosphorylated ERK.


HGF and IL-6 elicit cell type–specific activation dynamics of MEK

To distinguish receptor-specific and cell type–specific properties of ERK phosphorylation, we first focused on the time-resolved characterization of the input signal, the active, doubly phosphorylated MEK (ppMEK). We stimulated PMHs and HaCaT A5 cells with HGF (100 ng/ml) or IL-6 (100 ng/ml) and measured ppMEK by quantitative immunoblotting (fig. S1) (34). To describe the HGF- or IL-6–induced dynamics of MEK phosphorylation, we developed an input model of ordinary differential equations (ODEs) (Fig. 1, A and B, text S1, and table S1, A to I). We performed data-based model calibration for each cell type on the basis of χ2 statistics to evaluate how well the input MEK models captured the experimentally observed dynamics of ppMEK. We estimated the kinetic parameters of the mathematical model from our experimental observations.

Fig. 1 HGF and IL-6 elicit distinct response dynamics of ppMEK in PMHs and HaCaT A5 cells.

(A and B) Scheme of the input models in PMHs and HaCaT A5 cells. Arrows indicate biochemical reactions; round-headed arrows indicate activating reactions. Slashed circles represent degradation. (C and D) Calibration results of the PMH and HaCaT A5 models. Experimental data were acquired by quantitative immunoblotting (means ± SD, n = 3 biological replicates). Model trajectories are depicted in solid lines. The control represents no stimulation. See table S1 for model details.

In PMHs, HGF induced an increase in ppMEK, which peaked at about 10 min after HGF stimulation and then slowly decreased for up to 120 min (Fig. 1C). With IL-6, the ppMEK response was transient with a lower amplitude than that observed with HGF, which returned to basal levels by 120 min after stimulation. In HaCaT A5 cells, HGF induced a MEK activation response with a prolonged sustained appearance of ppMEK, whereas IL-6 triggered a transient, lower-amplitude response (Fig. 1D). The experimental data showed that HGF induced a greater amplitude of MEK activation compared with that induced by IL-6 and that this difference was independent of the cell type, suggesting that the ligands determined the amplitude of activation at this step of the MAPK pathway. However, HGF-stimulated MEK phosphorylation was sustained with possible oscillations in HaCaT A5 cells, indicating that signal duration and changes in other kinetic properties upon strong pathway activation may be cell type–specific.

Site-specific phosphorylation analysis reveals context-dependent activation of ERK

Because the dynamics of ppMEK activation displayed distinct features depending on the stimulus and the cell type analyzed, we explored the propagation of the signal to ERK. First, we determined the stoichiometric ratio of ERK1 to ERK2 by mass spectrometry (35) in combination with quantitative immunoblotting (34) (fig. S2, A and B) because previous studies have reported differences in isoform-specific signaling (3638). Indeed, in PMHs and HaCaT A5 cells, ERK2 exceeded ERK1 by factors of 2.0 and 2.6, respectively (fig. S2C). With this information, we converted the experimental data that represent fractions of total ERK1 or fractions of total ERK2, and the corresponding estimated errors, to relative fractions of the sum of total ERK1 and ERK2 by multiplying with the isoform ratio between ERK1 and ERK2 (fig. S2C). Thus, we show the abundances of phosphorylated forms as fractions of the sum of total ERK1 and total ERK2 and thereby maintain the relative ratio of all species to each other.

To elucidate the dynamics of site-specific ERK phosphorylation, we acquired densely sampled time-resolved data in PMHs and HaCaT A5 cells stimulated with HGF (100 ng/ml) or IL-6 (100 ng/ml) and determined the abundances of nonphosphorylated (np-), monophosphorylated (pT- or pY-), and double-phosphorylated (pTpY-) ERK1 and ERK2 by mass spectrometry with isotope-labeled peptide standards (35). We calculated data noise from the total degree of phosphorylation of ERK1/2 in each cell type by mathematical modeling (fig. S3, A and B, text S2, and table S2, A to G).

Despite their different abundances, the experimental data indicated nearly identical behavior of double-phosphorylated and activated ERK1 and ERK2 in each cell type (Fig. 2, A and B). The amplitude of HGF- and IL-6–mediated ERK1/2 activity (pTpY-ERK1/2) was qualitatively similar to that of the ppMEK signal induced by each ligand, with HGF triggering a stronger activation than IL-6. Similar to the profile of ppMEK, we observed a transient increase in pTpY-ERK1/2 upon stimulation with HGF and IL-6 in PMHs (Fig. 1C) and with IL-6 stimulation in HaCaT A5 cells (Fig. 1D, blue line), whereas HGF produced more sustained changes in pTpY-ERK1/2 in HaCaT A5 cells (Fig. 1D, green line).

Fig. 2 Data-based quantitative dynamic modeling results in distinct model structures of ERK phosphorylation in PMHs and HaCaT A5 cells.

(A and B) Top: Results of model calibration of the final PMH and HaCaT models. Solid lines depict trajectories of the calibrated model. Experimental data used for model calibration were acquired by mass spectrometry in biological quadruplicates, and the experimental error was estimated using a mathematical model based on ODEs (see table S2). The number of parameters and the goodness of fit of the model to the data are indicated. In PMHs, k16 and k19 have been estimated to the same value (see table S3C). Bottom: Scheme of the final ERK models in PMHs and HaCaT A5 cells using ppMEK as the input function. In the HaCaT model, two pTpY-ERK1/2–induced feedback loops, mediated by D and aPF, were included. Arrows indicate biochemical reactions; round-headed arrows indicate activating reactions. Slashed circles represent degradation. Dashed line indicates a reaction with delay. The models include phosphorylation and dephosphorylation of double-phosphorylated (pTpY-ERK), monophosphorylated (pT-ERK, pY-ERK), and nonphosphorylated (np-ERK) ERK. Because no isoform-specific parameters were included, the models represent both ERK1 and ERK2 as “ERK.” See table S3 for model details.

The two cell types exhibited a notable difference in the monophosphorylated ERK1 and ERK2 species. Compared to HGF-stimulated PMHs, which contained only small amounts of the monophosphorylated ERK1 and ERK2 species, HGF-stimulated HaCaT A5 cells exhibited a relatively high amount of pT-ERK1 and pT-ERK2 and a low amount of pY-ERK1 and pY-ERK2 (note the difference in y-axis scales in Fig. 2). Thus, our results indicated context-dependent signal propagation from MEK to ERK with characteristic variances in the site-specific phosphorylation of ERK.

Data-based quantitative dynamic modeling indicates distributive and cyclic ERK phosphorylation and dephosphorylation

To investigate the isoform- and site-specific phosphorylation dynamics, we developed a more detailed mathematical model of ppMEK-induced ERK phosphorylation, which includes all forward and reverse reactions and which we refer to as the “comprehensive” model, and calibrated the model with PMH and HaCaT A5 cell data (fig. S4, A and B, and table S3, A to F).

The model contained non-, mono-, and double-phosphorylated ERK1/2 species in response to ppMEK and initially consisted of all possible phosphorylation and dephosphorylation reactions using mass action kinetics, taking into account both distributive and processive ERK phosphorylation and dephosphorylation mechanisms (table S3D) based on previous reports (9, 10, 14, 3941). All phosphorylation reactions were controlled by ppMEK, whereas the dephosphorylation reactions were catalyzed by phosphatases with time-independent activities. The activities of the phosphatases were therefore not explicitly modeled but were included in the reaction rates of the dephosphorylation reactions. We performed model calibration with the PMH and HaCaT A5 cell data using χ2 statistics to evaluate how good the comprehensive model captured the experimentally observed dynamics of ERK (fig. S4B). Additionally, we assessed the ability of the mathematical model to sufficiently describe the qualitative features of the experimental data.

ERK in the model represents both ERK1 and ERK2. Because the experimental data revealed similar quantitative behavior for ERK1 and ERK2, we assumed nonisoform-specific kinetic rate constants, which reduced the number of free parameters. Although the comprehensive model captured the ERK dynamics in PMHs (fig. S4B, left), in HaCaT A5 cells the model failed to describe several qualitative features, including the peak of pTpY-ERK2 at 15 min and the basal amount of pY-ERK1 and pY-ERK2 (fig. S4B, right). Moreover, we also subjected the comprehensive model calibrated on PMHs and HaCaT A5 cells to structural identifiability analysis (42), which investigates how accurately the model parameters can be determined given the current model structure. Most of the parameters were estimated within a narrow range (fig. S4C). However, several parameters, for example, the rate constant for processive phosphorylation (kEpp), demonstrated a broad parameter range over several orders of magnitude at very small parameter values (10−8 to 10−6), indicating that the comprehensive model parameters were not structurally identifiable for both cell types.

We tested simpler structures of the comprehensive model to obtain a final model structure that described the experimental data obtained in the respective cell type with a minimal number of kinetic parameters. We first tested a model structure that only included a processive phosphorylation mechanism (39) and a distributive dephosphorylation mechanism (fig. S5A) (41), as previously suggested. However, the processive phosphorylation and distributive dephosphorylation model did not describe the PMH or HaCaT A5 cell experimental data. This model systematically underestimated the abundance of pTpY-ERK1/2 in both cell types (fig. S5B) and overestimated the amount of monophosphorylated ERK1/2 in PMHs (fig. S5B, left). In HaCaT A5 cells, the processive phosphorylation and distributive dephosphorylation model did not display any dynamics for pY-ERK1/2 (fig. S5B, right). For these reasons and on the basis of the likelihood ratio test (P < 0.05, likelihood ratio test for nested models in favor of the comprehensive model), we rejected this model structure.

We also performed a data-based model reduction using the parameter variances and excluded those reactions with corresponding rate constants that were compatible with zero (fig. S4C and table S3C). This reduced model accurately captured the ppMEK-induced transient behavior of pTpY-ERK1/2 and the dynamics of monophosphorylated ERK1/2 in PMHs (Fig. 2A), and thus, we refer to it as the “final” PMH model. The result of the model calibration suggested that dual ERK phosphorylation in PMHs only occurred through pT-ERK1/2 as an intermediate, whereas the deactivation of pTpY-ERK1/2 occurred either in a processive manner—directly from pTpY-ERK1/2 to np-ERK1/2 without the accumulation of any monophosphorylated species and without requiring any additional regulators or in a distributive manner—through pY-ERK1/2 to np-ERK1/2, resulting in cyclic phosphorylation and dephosphorylation (Fig. 2A). Upon model reduction, the HaCaT model described the measured ERK data with a goodness of fit that was similar to that of the comprehensive model (fig. S6, A and B; P > 0.05, likelihood ratio test for nested models in favor of the reduced model). Thus, we concluded that this model reduction was valid, and the reduction resulted in parameter values estimated within a narrow range (fig. S6C).

To improve the fit of the HaCaT model, we introduced two interconnected feedback loops emanating from pTpY-ERK1/2 (Fig. 2B and table S3, A to F). The first feedback loop consisted of the processive dephosphorylation of pTpY-ERK1/2 to np-ERK1/2 that required an additional dynamic regulator (D in the model in Fig. 2B), thus acting as negative feedback that enabled more precise modeling of the initial peak of pTpY-ERK1/2. Whereas the other phosphatases in the network operated in a time-independent manner, the activity of the dynamic regulator D is induced by pTpY-ERK. The second feedback loop catalyzed a processive activation of np-ERK1/2 to pTpY-ERK1/2, acting as positive feedback with an additional regulator (aPF in the model in Fig. 2B). Similar to the PMH model, the final HaCaT model structure suggested that cyclic phosphorylation and dephosphorylation could occur. Dual phosphorylation occurred through pY-ERK1/2 or through the processive reaction that required aPF, and dephosphorylation occurred through a species containing the phosphorylated threonine, either as pT-ERK1/2 or as pTpY-ERK1/2. Including these two feedback structures, this “final” HaCaT model accurately described the ppMEK-induced ERK phosphorylation dynamics (Fig. 2B). Using a likelihood ratio test for nested models, we verified that the model structure with feedbacks performs significantly better than the reduced model in describing the HaCaT data (P < 0.05 in favor of the final HaCaT model).

We performed model calibration using the data obtained in PMHs and HaCaT A5 cells for each of the final models (fig. S7, A and B). In these final models, all of the parameters were estimated within a narrow range (fig. S7, C and D, and table S4), indicating that each was a structurally identifiable model that was suitable for making predictions.

The network structure of the final HaCaT model suggested that ERK phosphorylation occurred through pY-ERK1/2 as an intermediate monophosphorylated form (Fig. 2B), which is in concordance with published results (43, 44). On the contrary, the reduced PMH model structure indicated that phosphorylation in PMHs occurred through pT-ERK1/2 as an intermediate monophosphorylated form (Fig. 2A) and suggested a high rate of conversion from pTpY-ERK1/2 to pY-ERK1/2 (table S4). To examine whether phosphorylation could also occur through the pY-ERK1/2 intermediate in PMHs, we tested two additional model structures that contained a distributive ERK phosphorylation mechanism through the pY-ERK1/2 intermediate and contained a positive or a negative feedback. We included dynamic feedback regulation affecting the dephosphorylation rate of pTpY-ERK1/2 to pY-ERK1/2 and tested whether this addition improved the fit to the experimentally observed amounts of pY-ERK1/2 at late time points at which ppMEK was already decreasing in PMHs (fig. S8, A and B). However, neither model structures sufficiently represented the experimental data, especially not the sustained increase in pY-ERK1/2 in PMHs, and thus both models were rejected. The qualitative disagreements were also reflected by high χ2 values of 717 (fig. S8A) and 727 (fig. S8B) compared to the χ2 value of 612 for the final PMH model (Fig. 2A). Therefore, we concluded that the final PMH model represented our best PMH model for describing the ERK phosphorylation and dephosphorylation mechanisms relevant in PMHs. Additionally, we showed that the final HaCaT model failed to describe the dynamics of the PMH data (fig. S9A) and that the final PMH model failed to capture the dynamics of the HaCaT data (fig. S9B and text S1). A constitutive processive phosphorylation reaction was unnecessary for model calibration and reduction, suggesting that a distributive or regulated processive phosphorylation mechanism accounts for ERK phosphorylation in both experimental systems.

DUSP6 and paxillin are involved in dual feedback regulation of ERK phosphorylation in HaCaT A5 cells

Our models suggested the presence of cell type–specific feedback regulation in HaCaT A5 cells that was absent in PMHs and was mainly triggered by HGF and not by IL-6. The induction of these feedbacks involved distinct time scales. The activation curve of the negative feedback (Fig. 3A, left) predicted a delay of 20.4 ± 1.9 min, suggesting transcriptional induction of the feedback regulator. The model predicted that the positive feedback (Fig. 3A, right) was activated immediately, indicating activation or recruitment of components that were already available.

Fig. 3 DUSP6 and paxillin are involved in feedback regulation in HaCaT A5 cells.

(A) HaCaT model predictions for the negative and positive feedback. Lines depict the simulated dynamics upon HGF and IL-6 stimulation and with no stimulation. (B) Experimental validation of DUSP6 as an induced negative feedback mediator in HaCaT A5 cells and not PMHs. Experimental data were acquired by qRT-PCR at the indicated time points (means ± SD, n = 3 biological replicates). Solid lines indicate model trajectories. (C) Experimental validation of paxillin as a nonenzymatic positive feedback mediator in HaCaT A5 cells. Cells were treated with scrambled siRNA or siRNA against paxillin for 48 hours and then stimulated with HGF (100 ng/ml). Paxillin, pTpY-ERK1/2, and total ERK1/2 were measured by quantitative immunoblotting. (D) Paxillin knockdown efficiency. Columns with error bars indicate means ± SD of all acquired data points (n = 15). (E) Model calibration of pTpY-ERK1/2 data. Experimental data (n = 3 biological replicates) were acquired by quantitative immunoblotting for treatment with scrambled siRNA or siRNA against paxillin. Solid lines indicate model trajectories. (F) Strength of positive feedback calculated as the kinetic rate constant of the positive feedback reaction estimated by model calibration.

DUSPs dephosphorylate ERKs, and some are transcriptionally induced by activated ERK within 30 to 60 min after stimulation with either cytokines or growth factors (4547). Thus, we screened for DUSPs that could participate in the negative feedback by measuring mRNA abundance using quantitative real-time PCR (qRT-PCR) (fig. S10A). We selected as candidates those DUSP transcripts that were specifically induced in HaCaT A5 cells upon HGF treatment. Transcripts encoding DUSP2, DUSP4, DUSP5, and DUSP6 were strongly induced in HaCaT A5 cells upon exposure to HGF. Because DUSP2 showed a transient induction, we rejected DUSP2 as a likely candidate. We also rejected DUSP5 because DUSP5 transcripts increased in both HaCaT A5 cells and PMHs exposed to HGF. We subjected DUSP4 and DUSP6 to model-based candidate selection (model extension is shown in table S5, A and B). Although both DUSP4 (fig. S10B) and DUSP6 (Fig. 3B) were induced with similar kinetics specifically in HaCaT A5 cells, the expression dynamics of DUSP6 exhibited the best fit with the model simulations (the model with DUSP6 had a χ2 value of 3.60, compared to the model with DUSP4 with a χ2 value of 23.79) (table S6); thus, we predicted that DUSP6 was the best candidate for this negative feedback regulator. Indeed, DUSP6 is localized in the cytoplasm and has high substrate specificity toward ERK compared to other MAPKs (26).

The final HaCaT model predicted a fast operating positive feedback loop that could be established by a scaffold. To identify potential candidates among the reported scaffolds of the MAPK pathway (48), we determined which scaffolding-encoding genes were expressed in each cell type (fig. S10C) and then attempted to knock down the scaffolds that were present in HaCaT A5 cells (fig. S10D). Of the scaffold-encoding transcripts present in HaCaT A5 cells, only paxillin and IQGAP1 were significantly reduced in abundance by knockdown (fig. S10D), indicating that the other scaffolds could not be effectively targeted or that they were constitutively expressed at a rate that could not be suppressed. We selected paxillin for analysis as a positive feedback regulator because paxillin has previously been reported as a scaffold protein involved in the HGF-induced MAPK pathway (25).

To investigate the importance of paxillin, we performed knockdown experiments and measured the amounts of paxillin, pTpY-ERK1/2, and total ERK1/2 by quantitative immunoblotting (Fig. 3C). The knockdown efficiency for paxillin was ~50% (Fig. 3D). By applying data-based model calibration, our model described the experimentally determined impact of paxillin knockdown on pTpY-ERK1/2 (Fig. 3E and table S7, A and B). Moreover, our model estimated a decrease in the strength of the positive feedback to 50% (Fig. 3F), which was consistent with the experimentally observed paxillin knockdown efficiency. Hence, we concluded that paxillin functions as a positive feedback regulator of ERK phosphorylation in HaCaT A5 cells.

Model-guided experimental design confirms the dynamics of DUSP6 and ERK after MEK inhibition

Having established identifiable quantitative dynamic models for the MEK/ERK module in PMHs and HaCaT A5 cells, we challenged the predictive power (49) of our models by comparing the model output to experiments with the MEK inhibitor U0126. To predict the most informative experiment, we used the final HaCaT model to simulate the DUSP6 expression after U0126 treatment at different time points (Fig. 4A and table S8, A to C). Because we found that DUSP6 mRNA reached its sustained plateau at 60 min after HGF stimulation (Fig. 3B), we exposed HaCaT A5 cells to U0126 60 min after HGF stimulation and compared the response to that predicted by the model. The model simulations predicted a decrease of DUSP6 mRNA to its basal amount 120 min after HGF stimulation if MEK was inhibited at 60 min. Congruently, whereas DUSP6 mRNA remained high in cells exposed to the dimethyl sulfoxide (DMSO) control after HGF stimulation, MEK inhibition 60 min after HGF exposure reduced DUSP6 abundance to basal levels by 120 min (Fig. 4B).

Fig. 4 Model-based prediction of the effect of MEK inhibition in HaCaT A5 is confirmed by experimental validation.

(A) Model simulation of the induced negative feedback with different time points of treatment with the MEK inhibitor U0126 at the indicated time points. Solid and dashed lines indicate the predicted model trajectories. (B) Experimental validation of the induction of DUSP6 mRNA by qRT-PCR. HaCaT A5 cells were stimulated with HGF (100 ng/ml), and 20 μM U0126 was added at t = 60 min. Open circles with error bars depict means and SD of the experimental data (n = 3 biological replicates).

We also tested the predictive power of the final PMH model. We compared the PMH model output in response to inhibition of MEK at the peak of pTpY-ERK1/2, which occurred 10 min after HGF exposure (Fig. 2A), to the response of PMHs exposed to HGF and then after 10 min exposed to U0126 (Fig. 5A). Here, we analyzed the effect of MEK inhibition on the dynamics of site-specific ERK phosphorylation by mass spectrometry. MEK inhibition by U0126 resulted in complete dephosphorylation of ERK 6 min after the addition of U0126 (16 min after HGF stimulation) (Fig. 5B, open circles). Although we observed a steep decrease for pTpY-ERK1/2, the amounts of pT-ERK1/2 and pY-ERK1/2 declined more slowly. As we found for HGF and IL-6 stimulation without MEK inhibition (Fig. 2A), the qualitative behavior between the ERK isoforms was similar. Because the kinetics of the inhibitory effect of U0126 on ppMEK are unknown, we included the decrease in ppMEK activity as an exponential decay (text S1 and table S9, A to D). Taking all other parameters from our calibrated model, we predicted the site-specific phosphorylation dynamics (Fig. 5B, solid lines). Our model was in good agreement with the experimental data of all measured ERK1/2 species. Thus, the final PMH model and the final HaCaT model of the MEK/ERK module accurately predicted the dynamic behavior of key components of this module after input perturbation, which was confirmed by independent, model-guided experimental validation.

Fig. 5 Model-based prediction of the effect of MEK inhibition in PMHs is confirmed by experimental validation.

(A) Setup of the validation experiment. After 10 min of stimulation with HGF (100 ng/ml), 20 μM U0126 was added to PMHs, and ERK phosphorylation dynamics were quantified by mass spectrometry. (B) Experimental data validation of the effect of MEK inhibition predicted by the PMH model. Data are shown as means and SD of the data (n = 3 biological replicates).

A preference for the double-phosphorylated form of ERK accounts for the higher amplitude of the PMH response to stimulation

To investigate the impact of the observed cell type–specific differences in the network structure on signal processing by the MEK/ERK module, we characterized the distribution profiles of the site-specific phosphorylated forms of ERK under three different scenarios in which the amount of ppMEK was low (baseline; untreated control), moderate (IL-6–induced), or high (HGF-induced), which corresponds to the experimental data representing these conditions (Fig. 6, A and B). For this analysis, we used the mathematical models that provided an analytical solution (data files S1 to S4) for the steady-state abundances of each phosphorylated form of ERK1/2 as a function of the amount of ppMEK. The values used to produce the analytical solution are provided in data files S2 and S4, and the models in SBML format are provided in data files S1 and S3. These analytical representations of the nonphosphorylated and phosphorylated ERK species in PMHs and HaCaT A5 cells revealed the dependence of the distribution of the phosphorylated forms of ERK on the amount of MEK activity.

Fig. 6 Information processing through ERK is characterized by the distribution of the phosphorylated forms of ERK.

(A and B) Distribution of phosphorylated forms of ERK at the baseline (basal), low (IL-6), and high (HGF) ppMEK range of the steady states based on the final PMH and HaCaT models. (C) Tables of criteria compliance or violation in PMHs and HaCaT A5 cells when individual parameters in the model were increased (×10) or decreased (×0.1) as indicated. The compliance or violation of the starting models (Model) for each of the three criteria is also indicated. Criteria are indicated between the PMH and HaCaT tables. (D and E) Dependency of ERK phosphorylated forms on the amount of ppMEK as determined by steady-state dose-response analysis. Continuous profiles for the three phosphorylated forms and the nonphosphorylated ERK were calculated for each model. Vertical dashed lines indicate the ppMEK concentrations reached by maximum HGF and maximum IL-6 stimulation and at baseline.

In PMHs, the analytical model predicted that most ERK1/2 was nonphosphorylated at baseline, indicating that the system was inactive when ppMEK was low (Fig. 6A). IL-6, which corresponds to a moderate amount of ppMEK, shifted the balance toward phosphorylated ERK1/2, with a preference toward pTpY-ERK1/2 at 37%. HGF, which corresponds to a higher amount of ppMEK, produced a greater increase in pTpY-ERK1/2 (65% of the total ERK1/2), and monophosphorylated ERK1/2 species were low (Fig. 6A). In contrast, the preferred phosphorylated species of ERK1/2 in HaCaT A5 cells upon stimulation with either HGF or IL-6 was pT-ERK1/2 (Fig. 6B). The fraction of pT-ERK1/2 exceeded pTpY-ERK1/2; thus, the model predicted that HaCaT A5 cells did not exhibit the bias toward pTpY-ERK1/2 that was predicted for PMHs.

There are, in principle, four possible mechanisms through which pT-ERK1/2 could accumulate: (i) increased phosphorylation from np-ERK1/2 to pT-ERK1/2, (ii) decreased dephosphorylation from pT-ERK1/2 to np-ERK1/2, (iii) decreased phosphorylation from pT-ERK1/2 to pTpY-ERK1/2, and (iv) increased dephosphorylation of pTpY-ERK1/2 to pT-ERK1/2.

To analyze at which point of the phosphorylation process the difference in the distribution profiles of the ERK1/2 phosphorylated forms arises, we systematically tested the effect of each kinetic parameter on the distribution of the phosphorylated ERK1/2 forms by varying the reaction rates by an order of magnitude (10-fold increase or 10-fold decrease) (fig. S11). We formulated three criteria for an unbiased characterization: “low basal phosphorylation,” representing the amount of np-ERK1/2 at the starting point; “dual phosphorylation preference,” representing a predominance of pTpY-ERK1/2 in response to MEK activation, and “efficient signal amplification,” representing the maximal amount of pTpY-ERK1/2 produced in response to MEK activation.

We first examined the PMH model and found that the best parameter set met all three criteria (Fig. 6C, model column). For PMHs, any parameter variation violated at least one criterion (Fig. 6C, orange panels), indicating that the original parameter set best fit the criteria. We therefore concluded that PMHs exhibited efficient signal amplification, which may arise through the combination of a low basal phosphorylation and a preference for the production of activated ERK1/2 (pTpY-ERK1/2) in response to stimulation.

In contrast to PMHs, the HaCaT model predicted that these cells failed to satisfy the dual phosphorylation preference criterion because of the relatively high amount of pT-ERK1/2 and low amount of pTpY-ERK1/2 (Fig. 6C). Thus, we individually varied the parameters to determine if any single change (increase or decrease by a factor of 10) enabled the system to meet all three criteria. An increase in the kinetic rate constant of pT-ERK1/2 dephosphorylation to np-ERK1/2 was the only parameter change that satisfied all three criteria (Fig. 6C and fig. S11B).

To investigate the effect of changes in the dephosphorylation rate of pT-ERK1/2, we performed model simulations of the dynamic behavior of ERK1/2 phosphorylation patterns (fig. S12). The simulations with a decreased dephosphorylation rate of pT-ERK1/2 showed an increased abundance of pT-ERK1/2 for the entire time course, a reduced amount of pTpY-ERK1/2, and an earlier peak of pTpY-ERK1/2 (fig. S12A). The dynamic model predicted that increasing the dephosphorylation rate of pT-ERK1/2 reduced the abundance of pT-ERK1/2 for the entire time course and only marginally affected the dynamics of pY-ERK1/2. This model simulation with the increased pT-ERK1/2 dephosphorylation rate also predicted an increased amount of pTpY-ERK1/2, which was in agreement with the distribution of the phosphorylated ERK forms predicted from the analytical model parameter variation analysis (fig. S11B). The peak timing of the pTpY-ERK1/2 dynamics was shifted to later time points. These results indicate that this particular kinetic parameter, which corresponds to reduced threonine phosphatase activity in HaCaT A5 cells, accounts for the high pT-ERK1/2 abundance observed in these cells when stimulated with HGF or IL-6.

Steady-state dose-response analysis reveals cell type–specific signal propagation through the MEK/ERK module

To characterize the cell type–specific distribution profiles of the phosphorylated ERK forms as a function of input strength (induced ppMEK activity), we performed a steady-state dose-response analysis as described previously by Suwanmajo and Krishnan (50). We used the analytical models to predict the change in each ERK species as a function of ppMEK. This analysis revealed a nonlinear relationship between ppMEK and ERK in PMHs (Fig. 6D) and HaCaT A5 cells (Fig. 6E). In both cell types, we observed three characteristic features: threshold, sensitivity, and saturation. In both of the PMH and HaCaT models, at baseline ppMEK concentration (black dashed line), the abundance of pTpY-ERK1/2 was low. In PMHs, the accumulation of pTpY-ERK1/2 showed a sigmoidal dependency on ppMEK concentration. In HaCaT A5 cells, a threshold on ppMEK concentration was more evident, whereas the amplitude of the increase in pTpY-ERK1/2 was less than that predicted for PMHs. At low ppMEK, corresponding to IL-6 stimulation (blue dashed line), the largest slope of the pTpY-ERK1/2 dose-response curve was observed in PMHs and HaCaT A5 cells, indicating a highly sensitive range. At a high ppMEK concentration, corresponding to HGF stimulation (green dashed line), pTpY-ERK1/2 was close to saturation for both models. For monophosphorylated ERK, the abundance of pY-ERK1/2 was low in both cell types, as was pT-ERK1/2 in PMHs. However, pT-ERK1/2 in HaCaT A5 cells showed a hyperbolic response with increasing ppMEK, comparable to pTpY-ERK1/2.

Specific ERK phosphorylation patterns exist in primary cells, cell lines, and hepatocellular carcinoma cells

On the basis of the analysis of our mathematical models, we postulated three criteria describing context-specific ERK phosphorylation patterns: basal phosphorylation status, dual phosphorylation preference, and efficient signal amplification. To investigate whether these patterns are tissue-specific or distinguish primary cells from immortalized cell lines, we analyzed the basal and HGF-induced ERK phosphorylation pattern in six cell systems by mass spectrometry.

We compared the ERK phosphorylation pattern in primary human keratinocytes and two HaCaT cell lines (Fig. 7A). We isolated primary human keratinocytes from breast tissue of a healthy human donor and grew the cells in serum-free medium supplemented with a mixture of defined growth factors including insulin and transforming growth factor–α (TGF-α) as described (51). Primary human keratinocytes were exposed to HGF (100 ng/ml) for 10 min or left untreated, and the ERK phosphorylation pattern was determined by mass spectrometry. Compared to either of the HaCaT cell lines, the primary human keratinocytes displayed less nonphosphorylated ERK at the basal state, with a relatively high amount of activated ERK (25% dual-phosphorylated). We predicted that this high basal ERK activation was due to the growth factors that are required for the cultivation of primary keratinocytes (51). In the primary keratinocytes, HGF induced an increase in the amount of pTpY-ERK1/2 and a decrease in the amount of pT-ERK1/2 and nonphosphorylated ERK1/2, suggesting that HGF induced a response with the dual phosphorylation preference. We performed the same experiment in the transformed human keratinocyte cell line HaCaT SD and the malignant human keratinocyte cell line HaCaT RT3. These cells were growth factor–depleted for 24 hours and stimulated with HGF (100 ng/ml) for 10 min or left untreated. Despite growth factor depletion, we observed substantial amounts of monophosphorylated ERK1/2 and, in the malignant cell line at the basal state, a relatively high amount of double-phosphorylated ERK1/2. This might be caused by autocrine secretion of growth factors by HaCaT cell lines as previously reported (31). Stimulation with HGF increased pT-ERK1/2 in HaCaT SD (P < 0.01, unpaired two-sided Student’s t test) and increased pTpY-ERK1/2 in both HaCaT SD and HaCaT RT3.

Fig. 7 Cell lines, primary cells, and tumor samples exhibit different distributions of the phosphorylated forms of ERK1/2.

(A and B) Experimentally determined distributions of the phosphorylated forms of ERK in primary cells and cell lines of the (A) skin and (B) liver. The abundance of the indicated forms of ERK was quantified by mass spectrometry in unstimulated cells and in cells stimulated with HGF for 10 min. Dots represent experimental data; bars indicate the average of the measurements (n = 2 to 3 biological replicates). Each cell type had specific culture requirements (see Materials and Methods). (C) Histology images (magnification, ×100) of tumor-free (left) and tumor-containing (right) liver samples from patients with cirrhosis who developed HCC. See fig. S13 for histology from other donors. (D) Site-specific ERK phosphorylation analysis by mass spectrometry (n = 9 patients). Quantified abundance of the indicated phosphorylated forms of ERK1/2 is shown as dots; the median is indicated by the horizontal line. Significance was determined by Student’s t test for paired samples: *P < 0.05.

To analyze the ERK phosphorylation pattern in hepatocytes (Fig. 7B), we isolated and cultivated (52) primary human hepatocytes from tumor-free liver tissue of two human donors diagnosed with colorectal liver metastases without additional pathological findings in the liver. In contrast to the keratinocytes, the unstimulated samples had very low amounts of mono- or double-phosphorylated ERK1/2, as expected due to the nontransformed status of the cells and the absence of growth factors in the cultivation medium. Stimulation of the primary human hepatocytes with HGF (100 ng/ml) for 10 min resulted in a strong increase in pTpY-ERK1/2 accompanied by only slight changes in monophosphorylated ERK species. This is consistent with the dual phosphorylation preference that we observed in PMHs. We performed the same experiment in the human hepatocellular carcinoma (HCC) cell line HepG2. Growth factor depletion without HGF stimulation resulted in low basal ERK phosphorylation. However, stimulation with HGF induced similar increases in pTpY-ERK1/2 and pY-ERK1 and a smaller increase in pT-ERK1, suggesting that these cells lack the dual phosphorylation preference that was characteristic of primary cells. Differences in the phosphorylation patterns in the basal state may reflect the need for growth factors for some of the cells (primary keratinocytes) or the production of autocrine factors that stimulate the MEK/ERK module [especially in the HaCaT cell lines (31)].

Because the primary mouse and human hepatocytes had low amounts of monophosphorylated ERK1/2 and the HepG2 cells, which are a liver carcinoma cell line, had a high amount of monophosphorylated ERK1/2, we hypothesized that the amount of monophosphorylated ERK might be associated with altered information processing in the MEK/ERK module that contributes to aberrant cellular proliferation. We obtained liver tissue samples from nine patients with cirrhosis who had developed HCC (fig. S13). We processed and analyzed by mass spectrometry paired samples from tumor tissue (Fig. 7C, left) and corresponding nontumor liver tissue (Fig. 7C, right). In tumor samples, pTpY-ERK1/2 was slightly increased, but the data were highly variable and the difference was not significant (P > 0.05). For the monophosphorylated ERK species, we found a significantly increased amount of pT-ERK1/2 in tumor samples compared to nontumor samples (P < 0.05) and no significant differences in the amount of pY-ERK1/2 (P > 0.05). Thus, increased amounts of monophosphorylated ERK, which reflect different information flow through the MEK/ERK module, may be indicative of dysregulated proliferation.


By combining an in-depth quantitative analysis of the dynamics of site-specific phosphorylation with data-based quantitative dynamic modeling, we revealed the complex organization of the MEK/ERK module of the MAPK pathway. In our experimental data, we observed differences in the response of the MEK/ERK module in HaCaT A5 cells and PMHs, which resulted in low amounts of pY-ERK1/2 in HaCaT A5 cells or the more sustained dynamics of pTpY-ERK1/2 in HaCaT A5 cells compared to the activation profile in PMHs. The mathematical models were based on previously established mechanisms of ERK phosphorylation and dephosphorylation, enabling dissection of cell type– and ligand-specific contributions to the regulation of ERK phosphorylation. Our study elucidated the impact of distinct network structures on cell type–specific ERK phosphorylation dynamics and quantitatively addressed the balance between the activation of kinases and phosphatases in shaping the distribution profiles of the phosphorylated forms of ERK. Our findings suggest that increased amounts of monophosphorylated ERK may be associated with cell transformation, both in cultured transformed cells and in tumor samples from patients with HCC.

Previous mathematical models of the MAPK pathway apply large models to a rather sparse data set (19) and focus on the overall structure of the pathway; smaller models with dense experimental data set facilitate the derivation of identifiable model parameters that provide insights into quantitative mechanisms (38). Modeling approaches using smaller models have stressed the importance of processive multisite phosphorylation of ERK, which can be achieved by either scaffold proteins or molecular crowding (53), resulting in faster signaling kinetics (6, 54). Other studies have demonstrated that ERK activation occurs in a distributive manner (10, 38, 55), suggesting that ERK phosphorylation mechanisms may be cell type–specific. For a detailed investigation of multistep ERK phosphorylation, time-resolved and site-specific phosphorylation data are necessary. Here, we combined such detailed quantitative measurements with mathematical modeling analysis and identified two model structures that described a distributive ERK phosphorylation mechanism in PMHs and HaCaT A5 cells. The power of the mathematical model is that it provides insights into the underlying biochemical mechanisms that cannot be intuitively derived from experimental data, such as the abundance of specific phosphorylated forms of ERK. Our model analyses revealed a cyclic network structure with phosphorylation and dephosphorylation mechanisms taking cell-specific routes. Such a cyclic structure might represent a common regulatory component of many signaling pathways, for example, as has been reported for the activation and deactivation of rhodopsin (56) or for the Kai circadian oscillator in cyanobacteria (57).

Our findings suggested that the network structure in HaCaT A5 cells includes negative and positive feedback loops. Our experiments indicated that DUSP6 was a likely candidate for a transcriptionally induced negative feedback regulator that shapes the activation dynamics of ERK (40). Our data also indicated that the positive feedback regulator was paxillin. Paxillin has been studied in detail in mouse inner medullary collecting duct epithelial cells (mIMCD-3). In mIMCD-3, paxillin is constitutively associated with MEK, and upon HGF stimulation, paxillin becomes phosphorylated and forms a complex with ERK, which brings both molecules in close proximity (25) and thereby can facilitate processive phosphorylation. A distributive phosphorylation mechanism can be converted to a processive mechanism (39). Theoretically, a distributive dephosphorylation can become more or completely processive in the presence of appropriate cofactors (41) or molecular crowding (39). These concepts are consistent with our identification of an induced processive ERK phosphorylation reaction (np-ERK to pTpY-ERK) in the HaCaT model, which may counter the induced negative feedback and enable the more sustained response kinetics.

Different phosphatase abundances might give rise to distinct distribution profiles of the phosphorylated forms of ERK, particularly concerning monophosphorylated ERK (58). The PMH model predicted that the distribution of phosphorylated forms of ERK was dominated by nonphosphorylated and double-phosphorylated ERK, especially with HGF, which stimulated the greatest response. PMHs exhibited some accumulation of pY-ERK1/2 (17.30%) upon IL-6 stimulation [the sum of 5.69% pY-ERK1 and 11.61% pY-ERK2 (Fig. 2A)]. However, this amount is low compared to the amount of pTpY-ERK1/2 (12.11% pTpY-ERK1 + 24.73% pTpY-ERK2 = 36.84%). The dephosphorylation of pT-ERK1/2 to np-ERK1/2 accounted for the relatively high abundance of pT-ERK1/2 in HaCaT A5 cells. Thus, our data and modeling indicated that differences in phosphatase abundances affect information processing through a phosphorylation-based signaling network. Such an effect has been shown for the ERK phosphatase DUSP1, the abundance of which is reduced in HCC. Further, alteration of the phosphatase and tensin homolog (PTEN) in the phosphatidylinositol 3-kinase (PI3K)–Akt network has been implicated in hematopoietic malignancies (5962).

The proliferation of PMHs during liver regeneration is primarily orchestrated by HGF (27). In contrast, HGF supports both the proliferation and migration of HaCaT A5 cells (29, 63). In PMHs, we observed a strong peak of ERK activation, which would be consistent with an HGF-induced proliferative response that would contribute to the liver mass regeneration following injury after two to three rounds of replication (64). However, high and sustained ERK activation can trigger cell cycle arrest and cellular senescence; this process is disrupted in tumor cells, which buffer ERK activation to an intermediate range to prevent senescence but promote proliferation (65). For example, in oncogenic Ras–transformed cells, attenuated ERK signaling prevents senescence and thereby enables cellular transformation (66). Congruently, HaCaT A5 cells exhibited a potentially oscillating ERK activation profile and a persistent, and potentially oscillating, increase in monophosphorylated ERK. These properties could be reflective of the transformed nature of these cells.

Analysis of tumor-containing and tumor-free liver samples from HCC patients supported the hypothesis that altered distributions of the phosphorylated forms of ERK may correlate with uncontrolled proliferation. It is generally assumed that activation of the MAPK pathway correlates with uncontrolled proliferation. Compared with adjacent normal tissue, the abundance of pTpY-ERK1/2, the active form of ERK1/2, was not statistically increased in the HCC samples that we analyzed. Instead, we detected higher pT-ERK1/2 in the tumor samples than in the respective normal liver tissue. Because excessive activation of ERK1/2 can lead to cell cycle arrest by induction of cell cycle inhibitors (6769), tumor cells may maintain pTpY-ERK1/2 at intermediate ranges. We propose that the pool of pT-ERK1/2 in transformed cells serves as a readily available intermediate for the dual-phosphorylated active form to facilitate unrestricted cell proliferation. Thus, the amount of monophosphorylated ERK might indicate altered information processing in the MEK/ERK module that can lead to the tumorigenic state and might constitute a more reliable cancer biomarker than double-phosphorylated ERK. In conclusion, data-based mathematical modeling of mass spectrometry data revealed differences in the molecular, topological, and kinetic properties of the MEK/ERK module between cell types, cell lines, and primary cells. Analysis of patient-derived materials suggested that increased relative amounts of monophosphorylated ERK may be a more useful indicator of persistent proliferative signaling than increases in the amount of active (double-phosphorylated) ERK.


Cell culture, stimulation, and time course experiments

The spontaneously immortalized nontumorigenic human keratinocyte–derived cell line HaCaT SD and the Ha-Ras–transformed tumorigenic cell lines HaCaT A5 and HaCaT RT3 (70, 71) were provided by P. Boukamp (Deutsche Krebsforschungszentrum, Heidelberg, Germany). HaCaT SD, HaCaT A5, and HaCaT RT3 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) containing high glucose (4.5 g/liter) (Lonza) supplemented with 10% fetal calf serum (FCS; Gibco) and 1% penicillin/streptomycin (Invitrogen). The medium for HaCaT A5 and HaCaT RT3 cells included G418 (200 μg/ml; Sigma). Primary human epidermal keratinocytes isolated from a healthy donor were provided by P. Boukamp (Deutsche Krebsforschungszentrum, Heidelberg, Germany). The cells were cultivated in DermaLife K Serum Free Human Keratinocyte Medium (CellSystems), consisting of DermaLife Basal Medium and LifeFactors [6 mM l-glutamine, 0.4% Extract P, 1.0 μM epinephrine, rhTGF-α (0.5 ng/ml), hydrocortisone hemisuccinate (100 ng/ml), rhInsulin (5 μg/ml), apo-transferrin (5 μg/ml), 0.02 to 0.15 mM calcium chloride] that were applied according to the manufacturer’s instructions. Cell lines were authenticated using Multiplex Cell Authentication by Multiplexion as described (72). The SNP (single-nucleotide polymorphism) profiles matched known profiles or were unique. The purity of the HaCaT cell line was validated using the Multiplex cell Contamination Test by Multiplexion as described (73). No Mycoplasma, squirrel monkey retrovirus, or interspecies contamination was detected.

PMHs were isolated from 8- to 12-week-old C57BL/6N mice (Charles River) housed at the DKFZ animal facility under a constant light/dark cycle, maintained on a standard mouse diet, and allowed ad libitum access to food and water. All animal experiments were approved by the governmental review committee on animal care of the state Baden-Württemberg, Germany (reference number A24/10). The isolation procedure was performed as previously described (28, 52). After isolation, hepatocytes were seeded at confluence (2 × 106 cells per 6-cm dish or 5 × 106 cells per 10-cm dish) in full medium [phenol red–free Williams E medium (Biochrom) supplemented with 10% (v/v) fetal bovine serum (Life Technologies), 0.1 μM dexamethasone, insulin (10 μg/ml), 2 mM l-glutamine (Sigma-Aldrich), and 1% (v/v) penicillin/streptomycin 100× (Life Technologies)] onto collagen I–coated cell dishes (BD Biosciences).

HaCaT cells and PMHs were cultured at 37°C, 5% CO2 and 95% relative humidity. For activation analysis of MEK1/2 and ERK1/2, PMHs, after cell adhesion, were washed with phosphate-buffered saline (PBS; pH 7.2 to 7.8) (PAN Biotech) and cultivated in serum-free cultivation medium [phenol red–free Williams E medium supplemented with 0.1 μM dexamethasone, 2 mM l-glutamine, and 1% (v/v) penicillin/streptomycin 100×] for 16 hours. PMHs were subsequently washed with PBS (PAN Biotech) and cultivated in serum-free cultivation medium depleted of dexamethasone for 6 hours. Fully confluent HaCaT cells were cultivated for 24 hours in FCS-free DMEM before treatment. Because primary human keratinocytes require serum-free medium including specific factors, the medium was not changed after 24 hours, and the experiment was performed using serum-free cultivation medium.

Primary human hepatocytes were isolated from two male patients who were diagnosed with liver metastasis derived from colorectal cancer. Informed consent of the patients was obtained according to the ethical guidelines of Charité—Universitätsmedizin Berlin. Primary human hepatocytes were isolated by a two-step collagenase perfusion technique according to Pfeiffer et al. (74). Primary human hepatocytes contained in the gained cell suspension were shipped in ChillProtec Plus (Biochrom) on ice overnight to DKFZ Heidelberg.

Primary human hepatocytes were cultivated in full medium as described for primary mouse hepatocytes (28), with an adhesion time of 6 hours. The human HCC cell line HepG2 (ATCC HB-806) was cultivated in phenol red–free DMEM containing high glucose (4.5 g/liter) (Gibco) supplemented with 10% FCS (Gibco), 1% penicillin/streptomycin (Gibco), 1% l-glutamine (Gibco), and 1% sodium pyruvate (Gibco). Cells were cultivated in serum-free medium for 16 hours before stimulation.

For time course experiments, PMHs were stimulated with recombinant mouse HGF (100 ng/ml) and HaCaT cells were stimulated with recombinant human HGF (rhHGF; 100 ng/ml) (both R&D Systems) or recombinant human IL-6 (R&D Systems) for the indicated times. To explore the influence of MEK1/2 inhibition, PMHs and HaCaT cells were stimulated with HGF (100 ng/ml) for 10 min and subsequently treated with 20 μM MEK1/2 inhibitor (U0126, Cell Signaling Technology) for the indicated times. Untreated or DMSO-treated cells served as a time-matched control.

Cell lysis, immunoprecipitation, quantitative immunoblotting, and SDS-PAGE

At the indicated time points after stimulation, cells were lysed using total cell lysis buffer [20 mM tris (pH 7.4), 150 mM NaCl, 1 mM EDTA (pH 8.0; AppliChem), 1% (v/v) NP40 (Roche Applied Sciences), 0.1% (w/v) sodium deoxycholate, 1 mM Na3VO4, 5 mM NaF, 4-(2-aminoethyl)benzenesulfonyl fluoride hydrochloride (AEBSF; 0.1 mg/ml), aprotinin (1 μg/ml)]. For ERK1/2 immunoprecipitation, total cell lysates were mixed with 10% SDS, incubated for 5 min at 95°C for denaturation, and subsequently diluted with 1× total cell lysis buffer. Protein A–Sepharose (GE Healthcare) and ERK1 antibody [ERK1 (K-23), Santa Cruz Biotechnology], which recognizes both ERK1 and ERK2, were added to the cell lysates; immunoprecipitations were performed overnight at 4°C, and immunoprecipitates were loaded on a 10% SDS–polyacrylamide gel electrophoresis (SDS-PAGE) gel. Recombinant human ERK1/2 (rhERK1/2) was loaded to indicate the corresponding endogenous ERK1/2 bands (Cell Signaling Technology; 0.5 μg of rhERK1 and 2 μg of rhERK2). Gels were stained with SimplyBlue SafeStain (Invitrogen) according to the manufacturer’s instructions. Bands corresponding to ERK1/2 were excised and subjected to mass spectrometry.

For the analysis of MEK1/2 and ERK1/2 activation and knockdown of scaffold proteins, total cell lysates were used. Total protein concentration was measured using bicinchoninic acid assay (Pierce), and lysates were loaded on a 10% SDS-PAGE gel, in case of MEK in a randomized order to reduce correlated errors (34, 75), separated by electrophoresis and transferred onto a polyvinylidene difluoride membrane (Millipore). The primary antibodies used were specific for phospho-MEK1/2 (Ser217 or Ser221) and total MEK1/2 (#9121 and #9122, respectively; Cell Signaling Technology), diphospho-ERK1/2 (M9692, Sigma), and total ERK1/2 (#9102, Cell Signaling Technology), β-actin (A5441, Sigma-Aldrich), paxillin (sc-365379, Santa Cruz Biotechnology), β-arrestin 2 (#3857, New England Biolabs), IL-17 receptor D (IL-17RD; ab96294, Abcam), IQGAP1 (AM20036AF-N, Acris), MP1 (sc-9136, Santa Cruz Biotechnology), and WDR83 (PAB7991, Acris). Secondary horseradish peroxidase (HRP)–coupled antibodies (anti-rabbit HRP, anti-goat HRP, protein A HRP) were purchased from Amersham Biosciences and Dianova and used for chemiluminescence detection using enhanced chemiluminescence substrate (GE Healthcare). Charge-coupled device camera–based signal detection was performed using an ImageQuant LAS 4000 system (GE Healthcare) and quantified using ImageQuant (GE Healthcare).

Mass spectrometry analysis with isotope-labeled peptide standards

For accurate and site-specific quantification of pT and pY of ERK1/2, mass spectrometry in combination with a targeted and label-based approach was applied. Isotopically labeled peptide/phosphopeptide one-source ratio standards were produced and applied as described in detail previously (35, 58, 76). These peptide standards were analogous to the activation motif–containing tryptic ERK1 and ERK2 peptides and were mixed with the tryptic in-gel digests after 37°C overnight incubation. Gel cutting, reduction, and alkylation, as well as digestion and peptide extraction, were performed as described previously (58). After peptide extraction from gel pieces, samples were purified for peptides by applying the STAGE-Tip method (77). Analyte peptide/standard mixtures derived from time course experiments performed in PMHs and HaCaT A5 cells used for model calibration were measured by nanoUPLC (nanoAcquity, Waters) coupled to an LTQ-Orbitrap XL mass spectrometer (Thermo Scientific). All other samples were measured by EASY-nLC 1000 (Thermo Scientific) coupled to a Q Exactive Plus Hybrid Quadrupole-Orbitrap mass spectrometer (Thermo Scientific). For both LC-systems, a precolumn setup and acetonitrile-based gradients (0 to 40% in <1 hour) were applied. To ensure a high number of survey scans (survey resolution = 60.000/70.000), the data-dependent acquisition methods were set to Top3 for Orbitrap-XL and to Top10 for Q Exactive. The degrees of ERK1 and ERK2 phosphorylation were analyzed by manual peak integration of endogenous and stable isotope–labeled standard peptides using Thermo Xcalibur Qual Browser version 3.0.63. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium ( through the PRIDE partner repository (78) with the data set identifier PXD002460. All measured abundances of the phosphorylated forms are shown as fractions of the sum of total ERK1 and total ERK2; for example, the abundance of pTpY-ERK1 was calculated aspTpYERK1(%)=pTpYERK1pTpYERK1+pTERK1+pYERK1+npERK1+pTpYERK2+pTERK2+pYERK2+npERK2By displaying all isoform- and phosphorylation-specific ERK species as fractions of the sum of total ERK1 and total ERK2, we maintain the relative ratio of all species. In our mathematical models, we treat ERK1, ERK2, and the phosphorylated forms thereof as individual variables.

Estimation of stoichiometric ratio of ERK and scaling of ppMEK data

For the estimation of the stoichiometric ratio of ERK1/2 in PMHs and HaCaT A5 cells, pTpY-ERK1/2 was measured in quadruplicates 10 min after stimulation with rhHGF (100 ng/ml). Antibodies (M9692, Sigma-Aldrich) with the same affinity to both double-phosphorylated ERK1 and ERK2 isoforms targeting the same sequence (HTGFLpTEpYVAT), as specified by the manufacturer, were used for immunoblotting (fig. S2A). This enabled the calculation of the ratio between pTpY-ERK1 and pTpY-ERK2 from the immunoblotting data. In addition, we performed immunoprecipitations of the same samples with antibodies directed against total ERK1/2 and analyzed the isoform-specific ERK ratios with respect to each phosphorylation site by mass spectrometry. On the basis of both ratios, we calculated the stoichiometric ratio between ERK1 and ERK2, and the sum of total ERK1 and total ERK2 was set to 100%. The estimated ratio (mean ± SD, n = 4) was used for the scaling of the mass spectrometry data. For activation analyses of MEK1/2, the maximal abundance of ppMEK was set to 5% on the basis of a previously published study (14) and other data points were scaled accordingly.

RNA isolation, RT-PCR, and qRT-PCR

For activation analyses of DUSPs, PMHs and HaCaT A5 cells were cultivated, incubated in serum-free medium, and stimulated as described above. At the indicated time points, RNA extraction was performed using the RNeasy Mini Kit (Qiagen) according to the manufacturer’s instructions. To obtain complementary DNA (cDNA) from RNA, the high-capacity cDNA reverse transcription kit (Applied Biosystems) was used according to the manufacturer’s instructions. qRT-PCR analysis was performed using LightCycler480 (Roche Applied Science). Samples were prepared with reagents of the LightCycler480 Probes Master Kit from Roche Applied Science. Specific primers were obtained from Eurofins MWG, and fluorescent hybridization probes [Universal Probe Library (UPL)] were obtained from Roche. For normalization of data, the geometric mean of the housekeeper genes HPRT, Ube2r2, and TBP was used. All primers and corresponding UPL probe numbers are listed in table S10, A and B; gene names are provided in table S10C.

Transfection with siRNA

HaCaT A5 cells were transfected with small interfering RNA (siRNA) oligonucleotides targeting human paxillin (SMARTpool: ON-TARGETplus PXN siRNA, L-005163-00-0005, Thermo Scientific), human β-arrestin 2 (SMARTpool: ON-TARGETplus ARRB2 siRNA, L-007292-00-0005, Thermo Scientific), human IL-17RD (SMARTpool: ON-TARGETplus IL-17RD siRNA, L-007946-02-0005, Thermo Scientific), human IQ motif containing guanosine triphosphatase–activating protein 1 (SMARTpool: ON-TARGETplus IQGAP1 siRNA, L-004694-00-0005, Thermo Scientific), human WD repeat domain 83 (SMARTpool: ON-TARGETplus WDR83 siRNA, L-014860-00-0005, Thermo Scientific), human late endosomal/lysosomal adaptor, MAPK, and MTOR activator 3 (SMARTpool: ON-TARGETplus LAMTOR3 siRNA, L-003572-00-0005, Thermo Scientific), and nontargeting siRNA (scrambled, SMARTpool: ON-TARGETplus Non-targeting Pool, D-001810-10-20, Thermo Scientific) serving as a control, using Lipofectamine RNAiMAX (Life Technologies) according to the manufacturer’s instructions.

Preparation of human liver tumor tissue

Tissue samples from liver tumor patients and corresponding nontumor liver tissue, as well as donor medical history and demographics (Table 1) and tissue histology, were provided by Charité—University Medicine Berlin. Informed consent of the patients for the use of tissue for research purposes was obtained corresponding to the ethical guidelines of Charité—University Medicine Berlin. Each liver tissue piece was cut, placed in tubes, and subsequently shaken twice for 30 s in a Precellys 24 homogenizer (Peqlab). Homogenized liver tissue was lysed in total cell lysis buffer [1% NP40, 150 mM NaCl, 20 mM tris-HCl (pH 7.4), 10 mM NaF, 1 mM EDTA (pH 8.0), 2 mM ZnCl2 (pH 4.0), 1 mM MgCl2, 2 mM Na3VO4, 20% glycerol, aprotinin (2 μg/ml), and AEBSF (200 μg/ml)] for immunoprecipitation. Lysates were rotated for 30 min at 4°C, sonicated, and centrifuged for 10 min at 14,000 rpm at 4°C. Four milligrams of the supernatant was incubated with the ERK1 antibody K-23 (Santa Cruz Biotechnology) and protein A–Sepharose (GE Healthcare) and rotated overnight at 4°C. Immunoprecipitated proteins were separated by SDS-PAGE. Gels were washed three times with water for 5 min and then stained for 1 hour with SimplyBlue SafeStain (Invitrogen) according to the manufacturer’s instructions.

Table 1 Donor medical history and demographics.
View this table:

Statistical analysis

Each experiment was repeated at least three times. Unless otherwise noted, data are presented as means ± SD.

Estimation of experimental error

We assumed a binomial distribution of our mass spectrometry data between 0 and 100%, depending on the total degree of phosphorylation (TDP) of all observed ERK species. For PMHs and HaCaT A5 cells, TDP was calculated asTDPi=pTpYi+pTi+pYi2and the corresponding measurement error (ME) was calculated asMEERKi,j=εERKi,j×ERKi,j×(100ERKi,j)with i = 1, 2 for ERK1 and ERK2 and j = 1, 2, 3, 4 for np-, pY-, pT-, and pTpY-ERK. The error parameters ε were estimated by applying the framework as described previously (49). We developed four ODE systems (table S9, A to G) that were solved by using the CVODES solver (79), and the model parameters were determined by maximum likelihood estimation. The final measurement error further considered the error of the determination of the stoichiometric ratio by propagation of uncertainty.

Mathematical modeling and parameter estimation

Development of the ODE models and model simulations were performed using the MATLAB toolbox PottersWheel version 3.0.11 (80). Because the abundances of MEK and ERK1/2 are similar (14, 38), all reactions were implemented as linear mass-action kinetics. Kinetic parameters were selected to be isoform-independent for ERK1 and ERK2. Therefore, the observed nonlinearities in information processing depended only on the model structure. Details of the developed models are provided in text S1, tables S1 to S3, S5, and S7 to S9, and data files S1 to S4. Parameters were estimated in logarithmic parameter space, based on a deterministic trust-region algorithm (81). For each parameter estimation round, up to 200 iterations were performed with a χ2 tolerance of 10−6 and a fit parameter tolerance of 10−6. For the indicated number of parameter estimation rounds, starting values of free kinetic rate constants were generated with a quasi-random number generator between 10−8 and 105. Initial values for model variables were given by the analytical steady state given by the model. Integration was carried out with the Fortran RADAU5 integrator based on an implicit Runge-Kutta method (82) of order 5 with an absolute error tolerance of 10−6 and a relative error tolerance of 10−3. Upon parameter estimation, χ2 values were calculated asχ2(P)=i=1N(yiy(ti;P)σi2)2with yi being data point i with SD σi and y(ti; P) being the model value at time point i for parameter values P. Because mass spectrometry–derived data were calculated as absolute numbers, no scaling or offset parameter was implemented.

Delay reactions in the mathematical model

The linear chain trick method (83) was applied to model a delay. Delay reactions capture a black box behavior with a certain time lag. The mean time delay was estimated as τ¯=Nk, with N being the number of implemented delay steps (N = 5, table S3A) and k being the reaction rate of the delay steps (k22, table S3C) (80).

DUSP screen data analysis

The transcript abundance of the DUSPs was measured in biological triplicates, and the mean and SD were calculated. In PMHs, there were two transcript variants for DUSP16 (DUSP16TV1, DUSP16TV2) and DUSP22 (DUSP22TV1, DUSP22TV2). In HaCaT A5 cells, there were three transcript variants for DUSP10 (DUSP10TV1, DUSP10TV2, DUSP10TV3). For these, the mean and SD were calculated over all replicates of all transcript variants (n = 6 for DUSP16 and DUSP22; n = 9 for DUSP10).

Paxillin data analysis

Paxillin was measured in biological triplicates. The signal intensities on a single blot were scaled on the mean of all measured signal intensities of each blot. The scaled signal intensities were averaged, and the SD was calculated (n = 3). The experimental knockdown efficiency was calculated by the ratio between the mean signal intensity of all paxillin siRNA–treated samples and the mean signal intensity of all control siRNA–treated samples, and the error considered the propagation of uncertainty.

Steady-state dose-response analysis

We analytically solved the ODE system to obtain the steady-state equations. For calculating the dose dependency of ERK on ppMEK, we formulated the steady-state equations as a function of ppMEK and the estimated parameter set (table S3C). We fixed the kinetic parameters according to the best models so that the steady-state equations were only dependent on ppMEK, which we varied from 0 to 5. For the analysis of criteria compliance, we either increased (×10) or decreased (×0.1) only the indicated parameter for an order of magnitude from the original value.

The sum of ERK1 and ERK2 was calculated for the characterization of the distribution of the phosphorylated forms. The values were quantified at three different ppMEK doses according to the ppMEK model in PMHs and HaCaT A5 cells. The ppMEK concentrations were ppMEK(PMH, baseline) = ppMEK(HaCaT, baseline) = 0.386 at baseline, ppMEK(PMH, IL-6) = 2.051, ppMEK(HaCaT, IL-6) = 1.166 at maximum IL-6 stimulation, and ppMEK(PMH, HGF) = 4.468, ppMEK(HaCaT, HGF) = 3.497 at maximum HGF stimulation.

Criteria formulation

For the characterization of the distribution of the phosphorylated forms of ERK, we established three formula-based criteria. Our first criterion considered the low ERK activation at its baseline (low basal phosphorylation), and it was fulfilled ifpTpYbaseline*pTpYbaselinewith the amount of basal pTpY-ERK1/2 in our original model denoted as pTpYbaseline and the amount of basal pTpY-ERK1/2 upon parameter variation (*) as pTpYbaseline.*

The second condition took into account the theoretical, statistical distribution based on the total phosphorylation degree and showed if the distribution of monophosphorylated ERK1/2 was low, indicating a dual phosphorylation preference (58). This condition, termed “dual phosphorylation preference,” was fulfilled ifpTHGF*pTHGFrpTIL-6*pTIL-6rpYHGF*pYHGFrpYIL-6*pYIL-6rwith the statistical distribution (r) pTri and pYripTir=pYir=TDPi×(100TDPi)100with the total degree of phosphorylation upon parameter variation (TDPi)*TDPi*=pTpYi*+pTi*+pYi*2and with i = 1, 2 for IL-6 and HGF. The third criterion, termed “efficient signal amplification,” considered a sufficient degree of signal transmission through the ERK phosphorylation circuit. It was fulfilled ifpTpYIL-6*pTpYIL-6pTpYHGF*pTpYHGF


Text S1. Quantitative dynamic modeling of the MEK/ERK module.

Text S2. Estimation of data noise of mass spectrometry data by data-based quantitative dynamic modeling.

Fig. S1. Quantitative immunoblotting of MEK and one-dimensional SDS-PAGE of ERK.

Fig. S2. Stoichiometric ratio of ERK1 and ERK2.

Fig. S3. Model-based estimation of experimental error.

Fig. S4. Data-based model calibration with the comprehensive ERK model.

Fig. S5. Data-based model calibration with the processive ERK model.

Fig. S6. Data-based model calibration with the reduced HaCaT model without feedbacks.

Fig. S7. Data-based model calibration with the final ERK model.

Fig. S8. Data-based model calibration with a modified PMH model with a dynamic feedback allowing ERK phosphorylation via pY-ERK1/2.

Fig. S9. Validation of the identifiable model structures by quantitative dynamic modeling.

Fig. S10. Screening of potential candidates of the negative and positive feedback loops in HaCaT cells.

Fig. S11. Distribution profiles of the phosphorylated forms of ERK after parameter variance.

Fig. S12. Results of model simulation after changing the rate constant of pT-ERK1/2 dephosphorylation.

Fig. S13. Histology of patient samples.

Table S1. MEK model: variables, input variables, additional variables, parameters, reactions, start value assignments, additional parameters, observed variables, and ODEs.

Table S2. Error model: variables, dynamic parameters, error parameters, reactions, observed variables, error structure, and ODEs.

Table S3. ERK model: variables, additional variables, parameters, reactions, observed variables, and ODEs.

Table S4. Parameter identifiability.

Table S5. Negative feedback: additional observed variables and additional parameters.

Table S6. Goodness of fit of model calibration of DUSP4 and DUSP6 data.

Table S7. Positive feedback: additional observed variables and additional parameters.

Table S8. Model prediction HaCaT cells: input variables, parameters, and reactions.

Table S9. Model prediction PMH: variables, parameters, reactions, and ODEs.

Table S10. Primers used in qRT-PCR experiments in PMHs and HaCaT cells and gene names.

Data file S1 (.xml format). SBML file of the PMH model.

Data file S2 (.txt format). Start value assignment for the PMH model.

Data file S3 (.xml format). SBML file of the HaCaT model.

Data file S4 (.txt format). Start value assignment for the HaCaT model.


Acknowledgments: We thank P. Boukamp for providing us with HaCaT cells and primary human keratinocytes and M. M. Mueller for conducting the initial experiments. Special thanks to S. Bonefas for the preparation of PMHs. We thank R. Merkle, K. Robichon, S. Chakraborty, and M. Uhrig for their support during the validation experiments. We thank A. Raue for his excellent technical support for establishing the phosphorylation degree-based error model. We thank F. Pinna for his support during the preparation of human liver tumor samples and M. Kießig for the preparation of histological slices and hematoxylin and eosin staining. Funding: Supported by German Federal Ministry of Education and Research (BMBF) grants to N.I. and M. Schilling (e:Bio collaborative research project “Systems Biology of Erythropoietin” SBEpo, 0316182B), to M.E.B. and W.D.L. (CancerSys network “LungSysII” 0316042A), to L.A.D. and P.L. (Virtual Liver Network 0315745), and to D.D., G.D., and D.S. (Virtual Liver Network 0315741) and by the German Center for Lung Research (DZL) to S.D. and U.K. (82DZL00404). Author contributions: M. Schilling and U.K. designed the project. W.D.L., L.A.D., S.D., and M. Schilling supervised the project. L.A.D. and S.D. performed cell culture, quantitative immunoblotting, and experiments for mass spectrometry measurements. B.H., M. Stepath, and M.E.B. performed mass spectrometry analysis with isotope-labeled standards. N.I. and B.A.K. conducted the measurements of mRNA by qRT-PCR. L.A.D., B.H., and M. Schilling analyzed the stoichiometric ratio of ERK1/2. N.I. performed data analysis, quantitative dynamic modeling, and further model analyses. P.L., A.V., D.D., G.D., and D.S. prepared the human liver samples. N.I. and M. Schilling wrote and all authors approved the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium ( through the PRIDE partner repository (78) with the data set identifier PXD002460.
View Abstract

Navigate This Article