Although the ERK pathway has a central role in the response of cells to growth factors, its regulatory structure and dynamics are incompletely understood. To investigate ERK activation in real time, we expressed an ERK–GFP fusion protein in human mammary epithelial cells. On EGF stimulation, we observed sustained oscillations of the ERK–GFP fusion protein between the nucleus and cytoplasm with a periodicity of ∼15 min. The oscillations were persistent (>45 cycles), independent of cell cycle phase, and were highly dependent on cell density, essentially disappearing at confluency. Oscillations occurred even at ligand doses that elicited very low levels of ERK phosphorylation, and could be detected biochemically in both transfected and nontransfected cells. Mathematical modeling revealed that negative feedback from phosphorylated ERK to the cascade input was necessary to match the robustness of the oscillation characteristics observed over a broad range of ligand concentrations. Our characterization of single‐cell ERK dynamics provides a quantitative foundation for understanding the regulatory structure of this signaling cascade.
Our increasingly detailed knowledge of cell signaling dynamics is enabling us to construct well‐constrained mathematical models that can illuminate regulatory mechanisms and lead to the in silico prediction of cellular responses (Wiley et al, 2003; Di Ventura et al, 2006). The wealth of quantitative data on the ERK activation pathway, its central role in eukaryotic cell decision processes and the complexity of its regulatory network structure have led to several mathematical models for this pathway (for review, see Orton et al (2005)). Some of these models have predicted the existence of oscillations in the pathway. However, despite numerous biochemical and imaging‐based investigations into its dynamics, oscillatory behavior of the ERK pathway has never been experimentally observed. Potential oscillatory behavior of the ERK pathway is intriguing because it was recently demonstrated that the NF‐kB and p53–MDM2 signaling pathways can display oscillations under certain conditions (Lahav et al, 2004; Nelson et al, 2004).
As the localization of ERK has been reported to reflect its state of activation (Lenormand et al, 1993), we used a fusion protein of ERK and green fluorescent protein (ERK–GFP) to follow its activation in individual cells. To simplify the analysis of ERK translocation, we used monomeric red fluorescent protein fused to a nuclear‐localization signal (mRFPnuc) to selectively visualize the nucleus and to allow us to track cells during long‐term live‐cell experiments.
When we stimulated cells with epidermal growth factor (EGF), we observed a rapid translocation of the ERK–GFP fusion protein into the nucleus. Significantly, in many cells the ERK–GFP fusion protein oscillated between the nucleus and the cytoplasm with a periodicity of ∼15 min (Figure 2C). The oscillations were asynchronous between neighboring cells and the mRFPnuc marker did not display any periodic oscillation. Interestingly, the frequency of the oscillations was essentially invariant over a wide range of conditions and EGF concentration. Oscillations required the continuous presence of EGF (Figure 2A) and were very sensitive to cell density, essentially disappearing at confluency.
To ensure that the oscillations were not an artifact of high expression levels of ERK–GFP fusion protein in our cells, we used flow cytometry to sort cells into different expression populations. We found that oscillation frequency was independent of the level of ERK–GFP fusion protein. By using ELISA and western blot analysis of phospho‐ERK levels in both parental and ERK–GFP‐expressing cells, we were also able to demonstrate oscillations in ERK phosphorylation that corresponded to the oscillations in ERK–GFP nucleo–cytoplasmic localization. Interestingly, overexpression of ERK–GFP suppressed the phosphorylation of the endogenous cellular ERK such that the net level of ERK phosphorylation (endogenous+transfected) remained constant, suggesting that the levels of ERK are not rate limiting for its activation in cells.
To understand the mechanistic basis of the oscillations, we constructed a mathematical model of the ERK pathway that included both the details of the biochemical reactions and transport processes between the cytoplasm and nucleus (Figure 5). Our initial parameters were based on literature values, but were adjusted to yield the observed oscillation characteristics. We found that a negative feedback loop between ERK and the input to the ERK cascade was essential to reproduce the waveform and invariant frequency aspects of the oscillations. Our model exhibited two Hopf bifurcations, and predicted that oscillations would occur over a range of two orders of magnitude of input strengths. Above and below a critical level of stimulation, however, the model predicted that oscillations would cease. These characteristics matched our experimental observations very closely.
Our initial model was based on the average oscillation properties of the cell population. To test its ability to predict the oscillation characteristics in a population of cells with cell‐to‐cell parameter variability, we assumed a log‐normal distribution of cellular parameters with a coefficient of variation 0.2 to be consistent with experimental measurements of the variability in protein expression between individual cells in a population (Spencer et al, 2009). We observed that model predictions for the population response to different input strengths were very similar to our experimental results and the model was able to quantitatively reproduce the robustness of the oscillations and the degree of cell‐to‐cell variability seen in our experiments. The model was also able to successfully predict the effect of inhibiting phosphatase activities on the fraction of oscillating cells as a function of EGF concentration.
Our study is, to the best of our knowledge, the first to experimentally demonstrate that ERK can oscillate between the nucleus and cytoplasm after cell activation. As illustrated by our analysis of ERK‐mediated negative feedback, ERK oscillatory dynamics can provide extremely useful information for constraining models of this important signaling pathway. We predict that oscillations could occur in other cell types depending on the inherent properties of their ERK cascades and on the existence of a negative feedback loop. We compared the parameters of our ERK cascade model with the parameters in other published models and found that our kinase and phosphatase abundances are within the range reported by others. However, in our model, the ratios of the enzymes’ Km values to their substrate concentrations are smaller compared with other published models.
Our investigation of the characteristics of the oscillatory behavior of the ERK pathway was greatly facilitated by using automated image analysis because of its ability to generate a very detailed ERK waveform. It could also follow hundreds of cells simultaneously at a time resolution of <1 min and provide data on associated behaviors, such as gene expression, cell migration and cell division. Live‐cell imaging of ERK dynamics should greatly facilitate the productive coupling of experiments to mathematical theory and help delineate the full regulatory structure of the MAPK cascade.
Activating the ERK cascade of human mammary epithelial cells with EGF resulted in a rapid and sustained oscillation of an ERK‐GFP fusion protein between the nucleus and cytoplasm.
Oscillations occurred even at ligand doses that elicited very low levels of ERK phosphorylation, were continuous throughout the cell cycle and were correlated with cycles of ERK phosphorylation detected biochemically in both transfected and non‐transfected cells.
Mathematical modeling revealed that negative feedback from phosphorylated ERK to the cascade input was necessary to match the robustness of the oscillation characteristics observed over a broad range of ligand concentrations.
When cell‐cell variations in levels of the different components of the ERK cascade was included in the mathematical model, we were able to accurately predict the population response to different levels of EGF stimulation and the effect of treating cells with phosphatase inhibitors.
Our increasingly detailed knowledge of cell signaling dynamics is enabling us to construct well‐constrained mathematical models that can illuminate regulatory mechanisms and lead to the in silico prediction of cellular responses (Wiley et al, 2003; Di Ventura et al, 2006). A good example is the ERK/MAPK pathway that exhibits complex and diverse dynamics and has been extensively studied over the past several decades (Seger and Krebs, 1995). Owing to its critical role in regulating cell proliferation, understanding ERK regulation is central to efforts to rationally design new antiproliferative drugs and other therapies (Roberts and Der, 2007). A potent regulator of ERK in epithelial cells is the epidermal growth factor (EGF) receptor that is activated by the binding of one of several cognate ligands (Citri and Yarden, 2006). Ligand binding induces activation of the intrinsic tyrosine kinase activity of the EGF receptor (EGFR), resulting in phosphorylation of both specific substrates and adaptor proteins, such as Shc, Grb2 and Sos. This then leads to the assembly of a membrane‐associated signaling complex that activates the ERK signaling cascade (Schoeberl et al, 2002). Phosphorylation of ERK by MEK results in its translocation to the nucleus where it phosphorylates transcriptional factors that drive the early response of cells to EGF (Lenormand et al, 1993).
The wealth of quantitative data on the MAPK activation pathway, its central role in eukaryotic cell decision processes, and the complexity of its regulatory network structure have led to several mathematical models for this pathway (for review, see Orton et al (2005)). These models have shown that the cascade structure and known regulatory features, such as positive and negative feedback loops and multisite enzyme phosphorylations, can contribute to diverse system properties, such as ultrasensitivity, bistability and oscillations. This response diversity could potentially enable the MAPK pathway to elicit distinct cellular responses depending on the cellular context. The MAPK activation cascade was shown experimentally to display strongly nonlinear or ultrasensitive responses more than a decade ago (Huang and Ferrell, 1996). The prediction of bistability has also been experimentally validated in studies that demonstrate its involvement in the generation of sustained responses to transient stimuli and irreversible cell‐fate decisions (Bhalla et al, 2002; Xiong and Ferrell, 2003).
Theoretical studies on the MAPK cascade have long predicted the existence of oscillations in the pathway. Negative feedback oscillations can be generated when the ultrasensitive MAPK cascade is coupled to a negative feedback loop wherein phosphorylated ERK inhibits the input to the cascade (Kholodenko, 2000). More recently, it has been shown that sequestration of kinases and phosphatases in the MAPK cascade by their substrates can potentially generate relaxation oscillations even in the absence of explicit feedback inhibition (Chickarmane et al, 2007; Qiao et al, 2007). However, despite numerous biochemical and imaging‐based investigations into its dynamics, oscillatory behavior of the ERK pathway has never been experimentally observed. Potential oscillatory behavior of the ERK pathway is intriguing because it was recently demonstrated that the NF‐κB and p53–MDM2 signaling pathways can display oscillations under certain conditions (Lahav et al, 2004; Nelson et al, 2004). Although the physiological significance of oscillations in these systems is not fully understood, it could be another regulatory layer in which biological information can be encoded.
Many cell types that are used in signal transduction research, such as HeLa cells, are used because they are easy to manipulate, even though they have limited physiological responses to growth factors such as EGF (Kinzel et al, 1990). In contrast, human mammary epithelial cells (HMECs) are similar to many epithelium‐derived cell types in that they require EGFR activation for proliferation and migration (Dong et al, 1999). The 184A1 strain of HMECs retains the EGFR‐dependent regulatory machinery of the primary cell type from which it was derived (Stampfer et al, 1993), and thus it is an excellent system for developing physiologically relevant models of the EGFR signaling pathway (Wiley et al, 2003; Hendriks et al, 2005; Shankaran et al, 2006). As ERK is a primary downstream effector pathway of the EGFR, we felt that an investigation of its potential oscillatory behavior in these cells was warranted. Here we show that the ERK cascade in HMECs displays very robust oscillatory behavior and that quantification of these oscillations can be used to help build accurate and predictive models of this important signaling pathway.
EGF induces oscillation of an ERK–GFP fusion protein between the nucleus and cytoplasm
After its activation by phosphorylation, ERK is rapidly translocated into the cell nucleus (Lenormand et al, 1993). Thus, we used the translocation of a fusion protein of ERK and green fluorescent protein (ERK–GFP) to assay ERK activation in individual cells. Cells were plated at a low density and brought to quiescence by overnight withdrawal of serum and EGF. After re‐addition of EGF, we observed a rapid increase of ERK–GFP fusion protein level in the cell nucleus (see Figure 1A). In many cells the ERK–GFP fusion protein oscillated between the nucleus and the cytoplasm with a periodicity of ∼15 min (Figure 1B and C and Supplementary Movies S1, S2 and S3). The oscillations were asynchronous between neighboring cells, and mRFP fused to a nuclear‐localization signal to serve as a nuclear marker did not display any periodic oscillation (Figure 1C).
The effect of EGF on nuclear translocation was specific to ERK and was correlated with its phosphorylation. Cells double‐labeled with ERK–GFP fusion protein and an mRFP marker without a nuclear‐localization signal showed relatively lower ERK–GFP fusion protein/mRFP ratio in the nucleus before EGF addition (Supplementary Figure S1). The addition of EGF caused a specific elevation of the ERK–GFP fusion protein/mRFP ratio in the nucleus and a simultaneous decrease in the ERK–GFP fusion protein/mRFP ratio in the cytoplasm with a periodicity of 12–15 min (Supplementary Figure S1). Immunofluorescence of fixed cells using antibodies against total ERK and phosphorylated ERK showed an increase in the levels of both species in the cell nucleus after EGF treatment, in agreement with previous studies (Supplementary Figure S2; Lenormand et al, 1993).
As shown in Figure 2A, oscillations required the continuous presence of EGF. After removal of EGF and blocking the EGFR with antagonistic mAb 225, cells underwent only a single additional nuclear translocation cycle. We then monitored the change in nuclear ERK–GFP fusion protein level over time in a field of cells grown at low density (Figure 2B). Each cell displayed a rapid oscillation pattern for the entire duration of the experiment (Supplementary Movie S4), and this pattern was sustained for more than 40 cycles (Figure 2C), at which point the experiment was terminated. There did not seem to be any relationship between the expression level of the ERK–GFP fusion protein and the oscillation pattern, and there was a low degree of synchrony between the cells in the population (Supplementary Movie S4). During this long experiment (>10 h), several cells underwent mitosis, allowing us to determine whether oscillations persisted throughout the cell cycle. Although there was a brief interruption in oscillations during mitosis itself when the nuclear membrane disappeared, on reformation of the nuclei, oscillations could once again be observed (Supplementary Figure S3). There was no obvious synchrony between the oscillations of the daughter cells. As we were following an asynchronous population of cells, continuous oscillations in all the cells suggest that ERK oscillations persist throughout the cell cycle.
Oscillations are independent of ERK–GFP expression levels
It has been suggested that the oscillations of NFκB, which have been observed after cell activation are a result of high expression levels of the GFP–NFκB fusion protein (Barken et al, 2005). To determine whether ERK oscillations could be due to high expression levels of our ERK–GFP construct, we used flow cytometry to sort out cells into ‘low’ and ‘high’ expression sets. The relative levels of ERK–GFP fusion protein versus endogenous ERK expression was then estimated by quantitative western blots (Figure 3A). This showed that the average level of ERK–GFP fusion protein in the population of low‐expressing cells was approximately two‐fold the endogenous level of ERK2. When the low‐expressing cell population was imaged after EGF addition, we observed oscillations regardless of the specific level of ERK–GFP fusion protein in individual cells (Figure 3B).
The translocation of ERK into the nucleus is postulated to be a consequence of ERK phosphorylation (Lenormand et al, 1993). Thus, we should also be able to detect oscillations in ERK phosphorylation states. We were unable to reliably image phosphorylated ERK in individual wild‐type cells over time using anti‐phosphoERK antibodies because of the low level of ERK expression (data not shown). However, we found that we could detect evidence of oscillating ERK phosphorylation in cell populations using either an ELISA (Figure 3C) or a western blot analysis (Figure 3E). The pattern of peaks and troughs, which we observed in the ERK–GFP‐expressing cells was very similar to that predicted from averaging their oscillating translocation pattern measured individually in multiple fields of cells (Figure 3D). We also measured the levels of phospho‐ERK in nontransfected HMECs after EGF treatment and observed a similar pattern (Figure 3C and E). Interestingly, we observed that the density at which the cells were grown had a pronounced effect on the population‐averaged oscillations of the ERK–GFP fusion protein (Figure 3D). The higher the cell density, the lower the fraction of oscillating cells. We confirmed the density dependency of ERK oscillations in nontransfected cells as well (Figure 3E and F). At densities of 1 × 104 cells per cm2, some degree of oscillation could be observed in the average cell population (Figure 3E), but at high cell densities, only an initial translocation could be detected (Figure 3F).
Quantitative analysis of ERK oscillation characteristics
To explore the apparent effect of context on ERK oscillations, we quantified the characteristics of the oscillations as a function of both cell density and EGF concentration. Using a waveform classification scheme (see Materials and methods and Supplementary information), we grouped oscillating cells into one of two groups: continuous oscillations (clean) or erratic oscillations (noisy). The difference between these two groups appears to be due to the relative magnitude of the oscillations, which allows them to be reliably detected by Fourier analysis (see Supplementary information). As shown in Figure 4A, cell density had a profound effect on the fraction of clean oscillators, which ranged from nearly 100% at very low densities to essentially none at confluence. The fraction of oscillating cells also decreased with EGF concentration, being maximal at 0.1 ng/ml and monotonically decreasing at higher concentrations (Figure 4B). The percentage of clean oscillators showed weak concentration dependence, decreasing from 54% at 0.1 ng/ml to 33% at 10 ng/ml. We found that the total level of ERK phosphorylation measured by an ELISA increased significantly with EGF concentration over the entire range used for the imaging experiments (Supplementary Figure S5). The ERK oscillations could be observed even at EGF concentrations that elicited only low levels of ERK phosphorylation (Supplementary Figure S5).
The shape and frequency of ERK oscillations, did not seem to be sensitive to cell density (Figure 4C, top panel) or EGF concentration (Figure 4D, top panel). The oscillation waveforms were distinctly asymmetric, with a rapid rise time of ∼6 min and a decay time of ∼9 min (Figure 4C and D, top panels). The oscillation time periods, determined using time domain analysis, varied between 11 and 23 min across the entire spectrum of EGF concentrations (Figure 4E). We estimated an average waveform periodicity of 15.9±4.3 min using Fourier Transform (FT) analysis and of 14.9±2.2 min using time domain analysis. Almost all individual oscillating cells displayed a decay time that was greater than the rise time irrespective of EGF dose (Figure 4F).
We also computed normalized oscillation amplitudes for the clean oscillators by dividing their absolute pulse amplitude with their nuclear ERK–GFP fusion protein level. The oscillation amplitudes of individual cells varied in a broad range under any given treatment (error bars in bottom panels of Figure 4C and D). The mean oscillation amplitude decreased significantly with cell density (Figure 4C) and showed a relatively weak proportional relationship with EGF dose (Figure 4D).
Mathematical model for nuclear–cytoplasmic ERK oscillations
To investigate the mechanism underlying ERK oscillations, we constructed a mathematical model for ERK activation and transport (Figure 5). Details of model construction are presented in Materials and methods section. Briefly, we chose the more general Huang–Ferrell formulation for the biochemical reactions in the ERK activation cascade (Huang and Ferrell, 1996; Qiao et al, 2007), in which each of the phosphorylation and dephosphorylation steps is modeled using explicit enzyme–substrate binding and catalysis reactions. The robustness of the ERK oscillation time period to changes in ligand concentration suggests that the underlying system is a negative feedback oscillator (Tsai et al, 2008). Hence, as in Kholodenko's treatment (Kholodenko, 2000), we included negative feedback from dually phosphorylated ERK to the cascade input (dotted red line in Figure 5), and chose our initial biochemical parameter values from that study.
In our preliminary analysis, we observed that our initial model formulation was in disagreement with one of the key features observed in our experiments. Our experimental data indicate that the oscillation time period is largely insensitive to the cellular ERK expression level (Figure 3B). To capture this feature, we had to make the extent of ERK phosphorylation insensitive to the ERK expression level. We achieved this by modeling the complex formation reactions between ERK and MEK species (reactions 7, 9 and 11; Figure 5) with rates that saturated with the ERK level becoming essentially independent of ERK concentration, rather than treating them as simple second order reactions. Biochemically, this could be accomplished by requiring scaffold proteins for ERK phosphorylation or simply making MEK rate limiting. We combined our biochemical model for the ERK cascade with the nuclear–cytoplasmic shuttling model given by Fujioka et al (2006) to model nuclear ERK dynamics (Figure 5). We modified biochemical rates and phosphatase abundances to match the nuclear ERK oscillation time period, rise and decay times measured in our experiments, while keeping the transport rates fixed. This constitutes our ‘base parameter set’ (Supplementary Table S1).
Oscillation characteristics predicted by the model for the base parameter set are presented in Figure 6. At an input strength E1tot=0.02 μM, oscillations of total nuclear ERK occur with a time period of 15 min, a rise time of 6 min and a decay time of 9 min. Oscillations in phospho‐ERK in the nucleus, the species relevant for transcription, are of an amplitude higher than that of total nuclear ERK (Figure 6A). The input strength E1tot can be interpreted as the total concentration of active MKKKK (Ras) available for stimulating the ERK cascade, and is expected to be proportional to the ligand concentration used (Tian et al, 2007). To examine the robustness of the oscillations to changes in the input, we determined the nature of the ERK response, and the maximum nuclear ERK‐PP levels for various values of E1tot (Figure 6B). We observed that the system exhibits two Hopf bifurcations, and that oscillations occur over a range of two orders of magnitude of input strengths from E1tot=2.5 × 10−3 to 2.7 × 10−1 μM. This large oscillatory regime is consistent with our experimental results for the effect of EGF concentrations on ERK oscillations. For E1tot <2.5 × 10−3 μM, the system stays in an ‘off’ state with negligible ERK activation. The lower Hopf bifurcation coincides with the threshold input strength beyond which ERK activation occurs. We observe a roughly linear dose response in the oscillatory regime (red dotted line, Figure 6B). For E1tot >2.7 × 10−1 μM, oscillations cease and we observe sustained ERK activation at a level that is relatively insensitive to the input strength. The graded dose response predicted for the oscillatory regime is consistent with our biochemical measurements of ERK phosphorylation under oscillatory conditions (Supplementary Figure S5). In the oscillatory regime the time period is predicted to vary from ∼7 to 19 min and the decay time is greater than the rise time for E1tot >∼8 × 10−3 μM (Figure 6C).
We also examined the effect of altering the total ERK expression levels in a six‐fold range on the oscillation characteristics at a fixed input strength E1tot=0.02 μM (Figure 6D). For constructing our base parameter set, we assumed that the total ERK concentration in wild‐type HMEC is 1.2 μM (expressed in terms of the cytoplasmic volume), a number that is based on estimates in other cell types (Huang and Ferrell, 1996; Brightman and Fell, 2000; Schoeberl et al, 2002; Fujioka et al, 2006). As the model is constructed as a means to explain the oscillatory behavior in ERK—GFP‐expressing cells, we set the ERK concentration in the base parameter set to 4.8 μM (Supplementary Table S1) four times higher than the assumed wild‐type level. When we varied the total ERK concentration from one to six times the assumed wild‐type level, oscillations were observed over the entire range of ERK expression levels tested. Although at the wild‐type ERK level (black line), oscillations occurred with a time period of ∼11 min, at two to six times the wild‐type level the oscillation time period was predicted to be ∼15 min. The robustness of the time period to ERK expression levels is consistent with our experimental observations (Figure 3B).
Qiao et al (2007) have recently shown that the MAPK cascade is capable of generating oscillations even in the absence of the explicit negative feedback loop. To examine whether a negative feedback loop is necessary to match our experimental observations, we removed the feedback and randomly sampled the biochemical parameter space of our model in a range of two orders of magnitude around the base parameter set. From our analysis of 10 000 random samples, we observed that ∼7% of the parameter sets exhibited Hopf bifurcations characteristic of oscillating systems. However, oscillations in the absence of feedback occurred over narrow ranges of the input strength, and with a time period that was extremely sensitive to the input. None of the oscillatory parameter sets matched our experimental observations in terms of the oscillation shape, time period robustness and width of the oscillatory regime (Supplementary Figure S6). This was true even when we sampled the parameter space in a range of two orders of magnitude around the parameter values listed in the study by Qiao et al (2007) (data not shown). In contrast, when we sampled the parameter space around our base set in the presence of negative feedback, we observed that 30% of the parameter sets were oscillatory, and were able to identify four parameter sets that matched the following stringent criteria: (i) oscillations over a range of two orders of magnitude of inputs, (ii) time period in the range of 5–25 min within the oscillatory regime, and (iii) decay time greater than the rise time for the majority of the inputs within the oscillatory regime. These additional parameter sets are listed in Supplementary Table S1. Overall, the explicit negative feedback loop was necessary to match our experimental observations. In the following sections, we present a more detailed analysis of our negative feedback model with the base parameter set. Qualitatively similar results would be obtained for any of the other four parameter sets listed in Supplementary Table S1.
Model predictions for oscillations in a cell population
We tested our model by evaluating its ability to predict the oscillation characteristics in a population of cells with cell‐to‐cell parameter variability. Our objective was to determine whether the model could accurately match our experimental results for the effect of EGF concentration on ERK oscillations (Figure 4B and D–F). We simulated cell‐to‐cell variability by sampling the inputs and parameters for each theoretical ‘cell’ from log‐normal distributions with a coefficient of variation (CV) of 0.2. This distribution is consistent with recently reported experimental measurements of the variability in protein expression between individual cells in a population (Spencer et al, 2009). Although the input strength E1tot will be proportional to the EGF concentration, the quantitative relationship between the two is not known. In our simulations, the total percentage of oscillating cells showed a biphasic dependence on the input strength with E1tot peaking at ∼0.01 μM (data not shown). Thus, we set the input strength of 0.01 μM to be equivalent to 0.1 ng/ml of EGF, as that was the concentration that induced maximal oscillations in our cells (Supplementary Figure S5). To estimate the fraction of clean oscillators, we defined a nuclear ERK oscillation amplitude cutoff of 0.07 μM—oscillations with amplitude lower than this value were deemed to be noisy. The amplitude cutoff was chosen so as to match the minimum oscillation time period determined in our experimental analysis of the clean oscillators (Figure 4E). In our simulations, the lowest observed time period was ∼6 min, but oscillations with time period lower than 10 min were of amplitude <0.07 μM.
Model predictions for the population response to different input strengths are in strong agreement with our experimental results (Figure 7). The total oscillating fraction and the percent of clean oscillators decreased when E1tot was increased from 0.01 to 0.2 μM (Figure 7A). As in the experiments, the oscillation time period, rise and decay times remained essentially unchanged with input strength (Figure 7B). The distribution of time periods was broader than in the experiments with the time period varying from 10 to 32 min (Figure 7C). It is possible that this is due to overestimation of the parameter variability in our simulations. Although protein abundances are expected to vary from cell to cell, it is unclear whether the rates of enzyme–substrate binding and catalysis would vary to the extent considered in our simulations. The model predicts an asymmetric oscillation shape as in the experiments with most of the cells displaying a decay time that is greater than the rise time (Figure 7D). Overall, the model is able to quantitatively reproduce the robustness of the oscillations and the degree of cell‐to‐cell variability seen in our experiments.
Model predictions for oscillation sensitivity to kinase and phosphatase levels
To examine the sensitivity of ERK oscillations to kinase and phosphatase levels, we constructed altered parameter sets by varying each of the three kinase and phosphatase abundances one at a time from one‐fifth to five times their value in the base parameter set (Supplementary Figure S7). We then quantified the fraction of oscillating cells in populations of 100 cells with parameters log‐normally distributed around each altered parameter set and inputs distributed around E1tot=0.02 μM. The fraction of oscillating cells showed the strongest dependence on the abundances of Raf and P1 (the MEK phosphatase) dropping to 0 when the Raf level was decreased to one‐fifth or the P1 level was increased to five times their respective base values. We observed modest sensitivity to the levels of MEK and P2 (the ERK phosphatase), and very little sensitivity to the levels of ERK and E2 (the Raf phosphatase) (Supplementary Figure S7).
To generate testable predictions, we simulated the effect of simultaneously inhibiting all of the phosphatases in the ERK cascade in a manner that would be expected when using an inhibitor such as sodium orthovanadate (Figure 8A). We simulated phosphatase inhibition by reducing all of the phosphatases in our model by the same fractional amount, and considered the population response to simulated input strengths of 0.02 and 0.001 μM. The former is expected to mimic the approximate effects of adding 1 ng/ml EGF. The input strength of 0.001 μM is 5% of the optimal value for the base set and is set to mimic the endogenous EGFR autocrine signaling that we observe in this cell type (Rodland et al, 2008). At an input strength of 0.02 μM, phosphatase inhibition is predicted to have a roughly monotonic inhibitory effect on ERK oscillations. Interestingly, at an input strength of 0.001 μM, ERK oscillations are predicted to show a biphasic dependence on phosphatase inhibition levels, increasing with phosphatase inhibition up to the 75% inhibition level and decreasing sharply thereafter.
We experimentally tested the model predictions for the effect of phosphatase inhibition on ERK oscillations by pretreating cells for 2 h with sodium orthovanadate before the addition of 0 or 1 ng/ml EGF (Figure 8B). We imaged ERK oscillations and quantified the fraction of clean oscillators as described in Materials and methods and Supplementary information. Untreated cells did not exhibit oscillations. However, pretreatment with 5 μM vanadate alone induced clean oscillations in approximately 30% of the cells. A higher concentration of vanadate (50 μM) resulted in a much smaller oscillating fraction. In experiments in which we added 1 ng/ml exogenous EGF, pretreatment with vanadate resulted in a dose‐dependent decrease in the oscillating fraction. These results are in qualitative agreement with our model predictions.
In our mathematical model, we assume that net ERK phosphorylation is not sensitive to increasing levels of ERK expression above a base value. This assumption was necessary to match the experimentally observed robustness of the oscillation time period, particularly its insensitivity to different levels of ERK–GFP expression. To test this assumption and to further validate our model, we treated both parental cells and ERK–GFP‐expressing cells with increasing concentrations of EGF and used western blot analysis to measure the levels of phosphorylation of both the endogenous ERK and the exogenous ERK–GFP fusion protein. As shown in Figure 8C and D, we observed that ERK phosphorylation had a similar EGF dose‐dependency in both cell types. Significantly, the total amount of phosphorylated ERK (ERK+ERK–GFP) in the transfected cells was similar to that observed in the parental cell despite the four‐fold level of ERK–GFP overexpression (Figure 3A). The suppression of endogenous ERK phosphorylation by ERK–GFP expression (Figure 8C) is probably due to substrate competition and indicates that cells have a limited capacity for ERK phosphorylation. This explains the relative insensitivity of ERK oscillations to the level of ERK–GFP expression, and serves as a further validation of our mathematical model.
Our study is, to the best of our knowledge, the first to experimentally demonstrate that ERK can oscillate between the nucleus and cytoplasm after cell activation. As illustrated by our analysis of ERK‐mediated negative feedback, ERK oscillatory dynamics can provide extremely useful information for constraining models of this important signaling pathway. In this study, we have used EGF activation of human epithelial cells as our model system. We have also observed that hepatocyte growth factor will induce oscillations in our cell type (our own unpublished observations) and have seen rapid oscillations of ERK–GFP fusion protein in mouse JB6 cells in response to fibroblast growth factor, but not EGF (T Weber and W Chrisler, unpublished observations). Interestingly, ERK oscillations in the JB6 cells occurred even at high cell densities. Thus, oscillations are likely to be a general property of the ERK pathway in many different cell types, although their characteristics will probably vary depending on the cell type and primary stimulus.
The observation that ERK oscillations have previously escaped detection might seem remarkable except for their characteristics. They are far more rapid than the sampling periods of most biochemical assays used to detect ERK activity. They are strongly inhibited by cell contact, at least in our cells, which would lower their chance of being detected using common protocols. The asynchronous nature of the oscillations would also tend to obscure their presence. Averaging nuclear ERK–GFP levels from even a small population of cells obscured oscillations beyond the first two cycles (Figure 3D). Finally, the use of transformed cells to investigate the ERK pathway might be a contributing factor because the dysregulation of signaling pathways that occurs during cell transformation could move the parameters of the ERK pathway out of the oscillatory regime.
Although ERK oscillations are remarkable for their persistence and regularity, whether they contain information that can cause differential cell responses is unclear. Extracellular signal‐regulated kinase is a potent activator of many nuclear transcription factors, and oscillations could be a means to selectively activate a subset of ERK‐responsive genes, analogous to oscillatory calcium signaling. In the case of calcium oscillations, information about stimulus dose can be encoded both in the amplitude and frequency of oscillations, which in turn have been proposed to control the level and specificity of gene expression (Dolmetsch et al, 1998). Unlike calcium oscillations, however, ERK oscillations do not display strong frequency or amplitude modulation in response to ligand dose. However, the strong dependence of the oscillation on cell density is consistent with it being a highly regulated process that could encode contextual information. It has been reported that different primary stimuli in PC12 cells can induce either transient or sustained activation of ERK and that these induce different cellular fates (Sasagawa et al, 2005; Santos et al, 2007). Conditions giving rise to oscillations are associated with an apparent sustained activation of ERK, whereas conditions that suppress oscillations give rise to transient ERK activation (Figure 3E and F). Thus, oscillation could be a mechanism underlying different cellular responses to persistent versus transient ERK activation. Although a direct role for ERK oscillations in controlling gene expression is intriguing, the oscillation could also simply be a consequence of the feedback control and the regulatory structure of the ERK pathway without directly encoding information. Experiments are underway to explore these different possibilities.
The existence of oscillatory behavior represents a valuable constraint that can be used to understand the regulatory structure of a biochemical system. Mathematical modeling becomes a particularly effective tool in this scenario. Our mathematical analysis of ERK activation and transport showed that an explicit negative feedback loop was necessary to match our experimental results. Although oscillations were possible in the absence of explicit feedback (Supplementary Figure S6), they only occurred over a narrow range of input strengths and had a periodicity that was strongly dependent on the strength of the input. For our base parameter set, the ERK cascade was extremely ultrasensitive (Hill coefficient, nH=393) in the absence of feedback. Introduction of negative feedback resulted in robust oscillations that occurred with a narrow range of time periods over a broad range of input strengths. This system is a pure negative feedback oscillator as shown in the study by Kholodenko (2000). We found 4 other parameter sets (Supplementary Table S1) that in the presence of negative feedback yielded similarly robust oscillations. Two of these involved ultrasensitive cascades (nH =5 and 17) that were rendered oscillatory due to the negative feedback (negative feedback oscillators). The remaining two solutions were more complex in that they exhibited relaxation oscillations as shown in the study by Qiao et al (2007) in the absence of feedback, but in a less than 1.3‐fold range of input strengths. The negative feedback loop was necessary to make these solutions robust oscillators, and hence observable under normal experimental conditions.
Whether oscillations occur in other cell types would depend on the inherent properties of their ERK cascades and the existence of the negative feedback loop modeled here. To further characterize the source of oscillations in our model and determine whether oscillations would be possible in other cell types, we compared the parameters of our ERK cascade model with the parameters in other published models for this reaction system (summarized in the Excel file in the Supplementary information). We observed that our kinase and phosphatase abundances are within the range reported in other models. However, for the negative feedback oscillators in our study (parameter sets P0, P3 and P4 in Supplementary Table S1) and in the study by Kholodenko (2000), the enzymatic reaction parameters are such that the kinases and phosphatases rapidly switch from an unsaturated state to a fully substrate‐saturated state during cascade activation and deactivation. Specifically, Km/S values—the ratios of the enzymes’ Km values to their substrate concentrations—are smaller in the oscillatory models compared with other published models. As noted in the study by Kholodenko (2000), this type of partial enzyme saturation contributes to strong cascade ultrasensitivity, and the possibility of oscillations in the presence of negative feedback. Although existing ERK models in the literature in general had higher Km/S values, we observed that most of these models were also ultrasensitive with Hill coefficients in the range 1.5–10. These ERK cascade models were constructed and parameterized to explain nonoscillatory, transient ERK activation in response to ligand stimulation. Hence, in these models the ERK cascade is couched within regulatory mechanisms that turn off the input to the cascade. It is possible that the nonlinearity of the ERK cascade is underestimated in this scenario, and that the true cascade parameters can only be determined in the absence of such regulatory mechanisms. Further, the negative feedback loop modeled in our study could also be present in other cell types, but only be evident when stronger regulatory mechanisms are absent. Thus, we suggest that oscillations may occur in many cell types, but might only be revealed using single‐cell observations under appropriate stimulus and cell culture conditions.
Our negative feedback model was able to correctly predict the effect of phosphatase inhibition on ERK oscillations (Figure 8A and B). Further, the model was able to quantitatively predict the oscillation characteristics in a population of cells in response to different ligand concentrations (Figure 7). We did not, however, attempt to model the effect of cell density on the oscillations, as the exact mechanisms by which density affects the ERK cascade remain to be established. At higher cell densities, we observed both a smaller ERK activation response and a transient response (Figure 3D). It is possible that intercellular interactions result in a sequestration or redistribution of the adaptor molecules required for ERK activation, thereby decreasing the input to the cascade. Additional regulatory mechanisms, such as an increase in the relative abundance of Dok to Shc in the model given by Sasagawa et al (2005), might be necessary to explain the transient nature of the ERK response at high cell densities.
Our investigation of the characteristics of the oscillatory behavior of the ERK pathway was greatly facilitated by our automated imaging method because it can generate a very detailed ERK waveform. It can also follow hundreds of cells simultaneously at a time resolution of <1 min and provide data on associated behaviors, such as gene expression, cell migration and cell division. Furthermore, the effects of introducing genes and siRNA into individual cells on ERK oscillations can be quantitatively evaluated and used to refine current models of ERK signaling. This should greatly facilitate the productive coupling of experiments to mathematical theory and help delineate the full regulatory structure of the MAPK cascade.
Materials and methods
Unless specified, all chemicals were purchased from Sigma. Antibodies against phosphoproteins were obtained from Cell Signaling Technology and secondary antibodies were from Jackson Immunoresearch Laboratories. Monomeric RFP was kindly provided by R Campbell (Department of Pharmacology, University of California at San Diego). Human mammary epithelial cell (HMEC) line 184A1 was obtained from M Stampfer (Lawrence Berkeley National Laboratory). The HMECs were maintained in DFCI‐1 medium (Band and Sager, 1989) as previously described (Wiley et al, 1998). Genes were introduced into these cells by retroviral transduction using virus from the medium of transiently transfected Phoenix cells (Kitamura et al, 1995).
Monitoring ERK translocation
The ERK1–GFP and mRFPnuc were subcloned into the retroviral vectors pBM‐IRESpuro (Garton et al, 2002) and pBM‐IRESblasticidin (the pBM‐IRES vector with blasticidin resistance substituted for puromycin resistance), respectively. The ERK1 was fused to GFP at its amino terminus, whereas monomeric RFP (Campbell et al, 2002) was fused to three copies of the nuclear‐localization sequence of Simian Virus 40 large T‐antigen at its carboxyl terminus by PCR using the EYFP‐Nuc vector as a template (Clontech). Cells expressing both markers were selected by antibiotic resistance and were found to grow at a rate indistinguishable from the parent cells.
Double‐labeled cells were seeded at varying densities onto a 24‐well Sensoplate coverlip (Greiner Bio‐One) and allowed to attach for 4 h at 37°C under 5% CO2. After 4 h, cells were washed and then incubated overnight with DFCI‐1 lacking EGF. Before imaging, the medium was replaced with 2 ml DFHB lacking EGF and bicarbonate, but containing 20 mM HEPES (pH 7.4). The cells were then placed onto a Nikon Eclipse TE300 inverted microscope equipped with an environment chamber (In Vivio Scientific, LLC) maintained at 37°C and allowed to equilibrate 1 h before imaging. Three random fields in each well were selected for imaging of both GFP (488 nm Ex; 508 nm Em) and mRFP (579 nm Ex; 604 nm Em) every 60 s with a Nikon Plan APO × 20/0.75 objective using a Retiga 1300 cooled CCD camera (QImaging) controlled by Volocity Acquisition (Improvision) software. EGF (Peprotech) was added in a 20‐μl volume to achieve the desired final concentration.
Epifluorescence images of mRFP‐tagged nuclei (red channel) were analyzed using the tracking utility in the Volocity software to track the locations of individual cell nuclei during the time course of the experiment. Nuclei were identified in the red channel using the SD Intensity function (2–100 limits), a medium noise filter, separating touching objects (size 0.1 μm2) and a pixel count of >100. Cell nuclei tracks starting from time t=0 and containing at least 150 continuous frames (total tracking time ∼150 min) were used for subsequent analysis. The mean ERK–GFP (green channel) intensity was measured in the regions occupied by the selected cell nuclei in each of the time frames. These measurements resulted in continuous nuclear ERK fluorescence versus time curves for individual cells.
Analysis of ERK phosphorylation
Cells were seeded in 100‐mm dishes in DFCI‐1 medium and allowed to attach for 4 h. The medium was then replaced with DFCI‐1 lacking EGF and cells were incubated overnight at 37°C. Before the experiment, the medium was replaced with DFHB lacking EGF and placed in a 37°C air incubator for 3 h. EGF was then added for the specified time and the reaction was terminated by a quick rinse with ice‐cold PBS containing 2 mM sodium orthovanadate (Na3VO4). Approximately 750 μl PBS containing 2 mM Na3VO4 and protease inhibitor cocktail set III (Calbiochem) was added to each dish and placed on ice. The cells were scraped into 2‐ml microcentrifuge tubes and centrifuged at 1500 r.p.m. for 3 min at 4°C using a swinging‐bucket rotor. The supernatant was aspirated and 20 μl of MT‐G lysis buffer (20 mM HEPES (pH 8.0), 1% Triton X‐100, 10% glycerol, 150 mM NaCl, 2 mM Na3OV4, 1 mM phenylmethylsulfonyl fluoride (PMSF), 1% aprotinin and protease inhibitor cocktail) was added to each tube and incubated on ice for 15 min. The cell lysates were cleared by centrifuging at 14 000 r.p.m. for 10 min at 4°C.
For western blot analysis, approximately 20 μg total protein of each sample was loaded onto a 10% Bis–Tris gel (Invitrogen), transferred to PVDF membrane and probed for phospho‐ERK. The p42/44 ERK phosphorylation was detected by immunoblotting with a 1:2000 dilution of mouse monoclonal phospho‐specific p42/44 ERK antibody with HRP‐conjugated goat anti‐mouse secondary antibody (Jackson Immunoresearch Laboratories) at 1:5000 dilution. Total p42/44 ERK was detected by immunoblotting with a 1:1000 dilution of rabbit polyclonal p42/44 antibody with HRP‐conjugated goat anti‐rabbit secondary antibody at a 1:5000 dilution. Quantification of p42/44 ERK phosphorylation and total ERK was performed using a Boehringer–Mannheim LumiImager and associated image analysis software.
For ELISA analysis, the phospho‐ERK1 (T202/Y204) sandwich ELISA kit from R&D Systems was used in a 96‐well format with slight modifications. Briefly, after cell lysis in MT‐G lysis buffer, approximately 40 μg of total protein was diluted to a 1:1 ratio with a lysis buffer supplemented with 8 M urea. The samples were vortexed and incubated at room temperature for 1 h and then centrifuged at 2000 × g for 5 min, and the supernatant was collected in a fresh tube. The samples were diluted to a final urea concentration of 1 M. Protein concentrations were measured in each sample using the BCA assay kit (Pierce). The ELISA assays were typically performed in triplicate for each sample and the results corrected for protein concentration, using 40 μg as the per‐sample standard.
Analysis of ERK oscillation characteristics
We used FT to estimate the frequency component of the oscillations associated with each cell. We also used a curve‐fitting approach to characterize the oscillatory waveforms in the time domain (see Supplementary information). We observed that most cells displayed a single oscillation frequency, although this could vary between different cells in a given field. By inspecting the waveforms in the time and frequency domain, we could classify cells as ‘clean’ oscillators, with a single, strong frequency component and a continuous train of at least 10 oscillations in the time domain (Supplementary Figure S4A); ‘noisy’ oscillators that display intermittent or very low amplitude oscillations, but with the FT spectrum showing a peak in the 10−3–2 × 10−3 Hz range (Supplementary Figure S4B); and nonoscillators with no dominant peak in the FT spectrum (Supplementary Figure S4C). Time domain analysis of the clean oscillators enabled us to quantify the time period, rise time, decay time and the oscillation amplitude (See Supplementary information).
Mathematical modeling of ERK oscillations
There are several existing models for the MAPK pathway, which vary widely with regards to model formulation, network structure and parameter sets (Orton et al, 2005). Further, two distinct mechanisms have been proposed for MAPK oscillations. The first prediction of oscillations in the MAPK cascade was made by Kholodenko (2000), who showed that cascade ultrasensitivity when combined with an explicit negative feedback loop from phosphorylated MAPK to the cascade input can lead to sustained oscillations. This model uses Michaelis–Menten rate expressions for the phosphorylation steps in the cascade. The negative feedback can arise either through the ERK‐mediated inhibitory phosphorylation of Sos (Langlois et al, 1995), or through the inactivation of Raf by ERK‐mediated hyperphosphorylation (Dougherty et al, 2005). Shvartsman and colleagues (Qiao et al, 2007) have recently shown that the MAPK cascade is capable of generating oscillations even in the absence of this explicit feedback, as kinase and phosphatase sequestration by their substrates can create implicit positive and negative feedback within the cascade. We constructed a mathematical model that would enable us to distinguish between these two possibilities.
For our analysis, we chose the more general Huang–Ferrell formulation for the ERK activation cascade (Huang and Ferrell, 1996; Qiao et al, 2007), in which each of the phosphorylation and dephosphorylation steps (reactions 1–10; Figure 6A) is modeled using explicit enzyme–substrate binding and catalysis reactions as opposed to Michaelis–Menten kinetics. This is the reaction system analyzed in the study by Qiao et al (2007), and is capable of oscillatory behavior in the absence of feedback. We modeled the complex formation reactions between ERK and MEK species (reactions 7, 9 and 11) with rates that saturated out with the ERK level, rather than treating them as simple second order reactions. We found that in the absence of this modification the ERK oscillation time period was extremely sensitive to the ERK expression level.
We set the total protein abundances for MKKK (Raf), MKK (MEK) and ERK for our model on the basis of previous estimates (Huang and Ferrell, 1996; Brightman and Fell, 2000; Schoeberl et al, 2002; Fujioka et al, 2006) and the assumption that ERK–GFP‐expressing cells have a four times higher ERK level than wild‐type cells. We included negative feedback from dually phosphorylated ERK to the cascade input (dotted red line in Figure 6A), and chose our remaining parameter values from the study by Kholodenko (Kholodenko, 2000) to form our initial parameter set for the ERK cascade. The negative feedback was modeled as an ERK‐PP‐induced reduction in the level of freely available active Ras (variable E1 in Figure 6A).
We combined our biochemical model for the ERK cascade with the nuclear–cytoplasmic shuttling model given by Fujioka et al (2006) to model nuclear ERK dynamics (Figure 6A). To make the biochemical and transport portions of the model compatible, we included the formation of complexes between unphosphorylated ERK and MEK (reaction 11, Figure 6A). The negative‐feedback model with the initial parameter set was capable of generating robust oscillations in nuclear ERK levels within specific ranges of the input strength (E1 in Figure 6A). We modified the biochemical rates and phosphatase abundances to match the oscillation time period, rise and decay times measured in our experiments, while keeping the transport rates fixed. This constitutes our ‘base parameter set’ (Supplementary Table S1).
Our negative‐feedback model consists of a single input, E1tot, which corresponds to the total level of externally imposed Ras activity. In addition to this, we have 42 biochemical parameters and 13 nuclear–cytoplasmic transport parameters, which together determine the dynamics of the 31 model variables. To examine the oscillation characteristics for the base parameter set (Figure 6), we chose a single value for the input strength E1tot and kept all of the parameters fixed as in the base set. We then varied either the input strength alone or the total ERK expression level as indicated in Figure 6. To identify other oscillatory parameter sets that matched our experimental observations, we randomly varied the biochemical parameters in a range of two orders of magnitude around the base parameter set, either in the presence or in the absence of the negative feedback loop. We fixed the MKKtot and ERKtot levels during this sampling, whereas MKKKtot was sampled randomly from one‐third to three times its value in the base parameter set. Finally, to examine the effect of cell‐to‐cell parameter variability on the oscillation characteristics, we sampled the input for each theoretical ‘cell’ from a log‐normal distribution with mean given by a specified value of E1tot and a CV of 0.2. The 42 biochemical parameters for each cell were sampled from a log‐normal distribution with mean given by the base parameter set and a CV of 0.2. Throughout our analysis, we only consider variations in the biochemical reactions of the cascade both for the sake of simplicity, and owing to the fact that the transport rates have experimentally determined values (Fujioka et al, 2006).
For any given set of biochemical parameters, we determined the existence of Hopf bifurcations by varying the input strength E1tot from 0 to 10 μM in the bifurcation analysis software AUTO (http://indy.cs.concordia.ca/auto/). AUTO reports the input strength values, if any, under which oscillations occur, and hence for a chosen input, we can immediately determine whether a cell would oscillate. To compute the oscillation profiles for a particular value of the input, we numerically integrated the model equations using MATLAB (Mathworks, Natick, MA) to determine the time period, rise and decay times. The total nuclear ERK concentration was used as the readout of the model to be compared with experimental observations. Detailed descriptions of the mathematical model and solution methodology are provided in the Supplementary information.
This study was funded by the Biomolecular Systems Initiative through the Laboratory Directed Research and Development Program at the Pacific Northwest National Laboratory (PNNL), a multiprogram national laboratory operated by Battelle for the US Department of Energy under Contract DE‐AC05‐76RL01830, and by NIH Grant 5R01GM072821‐03 to HR. A part of the research was performed using EMSL, a national scientific user facility sponsored by the DOE's Office of Biological and Environmental Research and located at PNNL.
Conflict of Interest
The authors declare that they have no conflict of interest.
Excel workbook containing a comprehensive comparison between the models and parameters presented in this paper and 7 previous models from the literature.
Video ‐ S1
Cells expressing ERK‐GFP prior to addition of EGF
Video ‐ S2
Cells expressing ERK‐GFP after addition of EGF
Video ‐ S3
Same sample as Movie S2, but imaging the red channel
Video ‐ S4
EGF‐induced ERK oscillations in a field of cells over a 2 h time period.
Supplementary figures S1–S7, supplementary methods and supplementary tables S1–S2
This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.
- Copyright © 2009 EMBO and Nature Publishing Group