Skip to main content
Download PDF

Regulation of ERK-MAPK signaling in human epidermis

BMC Systems Biology volume 9, Article number: 41 (2015) Cite this article

Abstract

Background

The skin is largely comprised of keratinocytes within the interfollicular epidermis. Over approximately two weeks these cells differentiate and traverse the thickness of the skin. The stage of differentiation is therefore reflected in the positions of cells within the tissue, providing a convenient axis along which to study the signaling events that occur in situ during keratinocyte terminal differentiation, over this extended two-week timescale. The canonical ERK-MAPK signaling cascade (Raf-1, MEK-1/2 and ERK-1/2) has been implicated in controlling diverse cellular behaviors, including proliferation and differentiation. While the molecular interactions involved in signal transduction through this cascade have been well characterized in cell culture experiments, our understanding of how this sequence of events unfolds to determine cell fate within a homeostatic tissue environment has not been fully characterized.

Methods

We measured the abundance of total and phosphorylated ERK-MAPK signaling proteins within interfollicular keratinocytes in transverse cross-sections of human epidermis using immunofluorescence microscopy. To investigate these data we developed a mathematical model of the signaling cascade using a normalized-Hill differential equation formalism.

Results

These data show coordinated variation in the abundance of phosphorylated ERK-MAPK components across the epidermis. Statistical analysis of these data shows that associations between phosphorylated ERK-MAPK components which correspond to canonical molecular interactions are dependent upon spatial position within the epidermis. The model demonstrates that the spatial profile of activation for ERK-MAPK signaling components across the epidermis may be maintained in a cell-autonomous fashion by an underlying spatial gradient in calcium signaling.

Conclusions

Our data demonstrate an extended phospho-protein profile of ERK-MAPK signaling cascade components across the epidermis in situ, and statistical associations in these data indicate canonical ERK-MAPK interactions underlie this spatial profile of ERK-MAPK activation. Using mathematical modelling we have demonstrated that spatially varying calcium signaling components across the epidermis may be sufficient to maintain the spatial profile of ERK-MAPK signaling cascade components in a cell-autonomous manner. These findings may have significant implications for the wide range of cancer drugs which therapeutically target ERK-MAPK signaling components.

Background

The epidermis is an epithelial tissue which forms the outermost layer of the skin (Fig. 1a) and performs an essential role in protecting an organism against environmental perturbations [1]. Interfollicular human epidermis consists of keratinocytes arranged in a spatial gradient of differentiation across the deep-to-superficial axis of the tissue (Fig. 1b). Epidermal barrier function is dependent upon changes to keratinocyte biochemistry and morphology that occur during differentiation [1, 2] and discrete tissue layers (Fig. 1b) are defined by histological features that reflect underlying molecular changes [1, 3–7]. A variety of signaling proteins are also modulated across the epidermis to coordinate keratinocyte terminal differentiation (Fig. 1c), including extracellular matrix ligands and associated plasma membrane receptors [8–10], and calcium (Ca2+) signaling components [11, 12]. Of particular interest for this study, increasing extracellular Ca2+ ion concentration is known to promote keratinocyte terminal differentiation in vitro [13, 14], and a Ca2+ gradient is maintained across the depth of the epidermis, increasing from the basal layer to the outermost granular layer (similar to ‘Superficial Signals’ in Fig. 1c) before decreasing across the transitional layer [11, 12].

Fig. 1
figure 1

ERK-MAPK Signaling Within Human Epidermis. a Epidermis is the outermost tissue layer of the skin with an essential role in protection from the environment. Epidermal barrier function is established and maintained by keratinocytes which undergo large biochemical and morphological changes during keratinocyte terminal differentiation. This establishes a spatially-regulated keratinocyte differentiation gradient across the depth of the epidermis, between hair follicles (within interfollicular epidermis). b Differentiating keratinocytes are pushed towards the superficial surface of the epidermis by proliferation within the basal layer. As this occurs, keratinocytes undergo terminal differentiation, establishing a spatiotemporal differentiation gradient across the depth of the epidermis. c The effect of tissue structure on paracrine/endocrine signals, and differentiation-associated changes in the abundance or activity of scaffold co-factors establish signal gradients across the depth of the epidermis. The gradient of Ca2+ within the epidermis is similar to the ‘superficial signals’ example; however, it peaks just prior to the transitional layer, rather than within the outermost superficial layers. d A simple representation of the canonical ERK-MAPK signaling cascade with: inputs to Raf-1 from extracellular calcium (Ca2+; activating) and plasma membrane calmodulin (CaM; inhibiting) modulated by cellular position along the keratinocyte differentiation gradient; the signal transduction cascade through Raf-1, MEK-1/2 and ERK-1/2; negative feedback from phospho-ERK-1/2 to phospho-Raf-1; and nuclear phospho-ERK-1/2 promoting its own dephosphorylation. Further details on these interactions are given within the Materials and methods. Nodes drawn in grey are not explicitly modeled, as they were not measured experimentally. A more comprehensive reaction kinetic scheme is given within Additional file 7: Figure S1 in reaction_network.png

Full size image

The canonical extracellular signal regulated kinase (ERK) cascade of the mitogen activated protein kinase (MAPK) family has been implicated in the regulation of keratinocyte differentiation in vivo [15–17] and in vitro [18, 19]. Signaling through ERK-MAPK integrates and mediates the effects of epidermal growth factor receptor (EGFR) [19, 20], integrin [21, 22] and calcium signaling [18, 23]. These input stimuli are localized to the plasma membrane where they regulate the conversion of Ras-GDP to active Ras-GTP. Generation of Ras-GTP promotes signaling through a sequential cascade of kinases which are activated by phosphorylation prior to phosphorylating their own downstream targets, progressing through activation of Raf-1 and B-Raf dimers, to MEK-1/2 and then ERK-1/2 (Fig. 1d). Phosphorylated ERK-1/2 then act upon a large number of proteins including several transcription factors to activate gene expression and initiate a cellular response [16, 24]. The sub-cellular localization of activated ERK-MAPK components is also important for determining specific cellular responses to the wide range of signals that influence ERK-MAPK signaling [25, 26].

The ERK-MAPK cascade has been extensively characterized using in vitro models [27–30] and computational methods [25, 31–33], typically over a time course of minutes-to-hours following addition of a mitogenic stimulus [27] or some other perturbation [29]. A number of studies have modeled intracellular signaling during developmental pattern formation, primarily using Drosophila and other malleable systems [34, 35]. Within human tissues however, the temporal and spatial dynamics of ERK-MAPK signaling, and the role that it may play in regulating cellular processes to help maintain tissue homeostasis are not well understood. Keratinocytes take approximately two weeks to traverse the epidermis as they undergo terminal differentiation [36–38], raising the question as to how ERK-MAPK signaling operates over this much longer timescale in situ. Furthermore, cells grown in vitro are exposed to the different environmental cues than cells in their in vivo tissue, and this can have significant effects on phenotypic behavior and intracellular signaling. For example, epidermal keratinocytes grown in vitro undergo abnormal differentiation expressing proteins more commonly associated with wound-healing and basal cell carcinoma [39–41].

In this study we therefore determined to characterize ERK-MAPK signaling in human epidermis in situ, by measuring abundance of signaling pathway components and mathematically modeling the spatiotemporal dynamics of the signaling cascade. We used immunofluorescence to measure the abundance of phospho-Raf-1, phospho-MEK-1/2, phospho-ERK-1/2 and calmodulin (CaM) within healthy human skin samples. A keratinocyte’s position within the epidermis, along the deep-to-superficial axis of the tissue, directly reflects its relative stage of differentiation. We used this to develop a computational framework for transforming the image data into a quantitative format, which allowed us to combine observations across multiple experiments. Our results show gradual, coordinated increases in the abundance of ERK-MAPK phospho-proteins over the depth of the skin, over much longer timescales than those observed for in vitro studies. Statistical analyses highlighted the importance of considering spatial position within the tissue when examining relationships between the measured variables and suggested canonical interactions are active during keratinocyte differentiation. Thus, we constructed a mathematical model of ERK-MAPK signaling interactions derived from the literature, using a normalized-Hill differential equation approach [42, 43], with which to analyze our data. The purpose of this model was to test whether Ca2+ signaling components which vary over the depth of the epidermis were sufficient to drive signaling through canonical ERK-MAPK interactions and reach different steady-state activation levels, in a manner which was consistent with the extended spatiotemporal phosphorylation profiles observed within our data. Our results suggest that canonical ERK-MAPK interactions are operating, as overall pathway activity is modulated through the depth of the epidermis.

Results

Spatial transformation of immunofluorescence data using histological landmarks reveals distinct spatial profiles of in situ ERK-MAPK phosphorylation within human skin

We used confocal microscopy imaging of immunofluorescence-labeled human epidermis to measure the abundance of ERK-MAPK phospho-proteins and calmodulin within interfollicular keratinocytes from samples of human skin (Additional file 1: Figure S2, Additional file 2: Figure S3, Additional file 3: Figure S4, Additional file 4: Figure S5 in pRaf1.png; pMEK.png; pERK.png; CALM.png). Regions of interest were manually selected within the image data using a GUI which displayed the images and samples at different z-positions and tracked input co-ordinates through ‘mouse clicks’. This provided a relatively quick approach to sample the image data across multiple patients and targets, while distinguishing distinct subcellular domains of the cytoplasm, nucleus (Fig. 2a), and where possible the plasma membrane (Additional file 4: Figure S5 in CALM.png).

Fig. 2
figure 2

Transforming fluorescence image data in to a spatially-conditioned, quantitative format for mathematical analysis. a Image data were obtained by immunofluorescence labeling and confocal microscopy, and different sub-cellular localizations (cytoplasm and nucleus) were manually sampled using a graphical user interface. The position and orientation of the cell selected for surface rendering (inset) has been highlighted (orange lines). b Epidermal tissue layers that could be distinguished using label-independent criteria were demarcated. The relative position of each sample was normalized within the layer using linear interpolation, then added to a whole integer which distinguished tissue layers (ilayer; Table 1), such that the normalized distance, dnorm = ilayer + (d1/(d1 + d2)). c The sampled fluorescence intensity data (grey dots) underwent loess smoothing (blue line). Spatial conditioning allowed data from Patients One (red), Two (green) and Three (blue) to be directly compared. d Protein abundance data from specific sub-cellular compartments were compared to a normalized-Hill differential equation model, with a literature derived network structure (described in Fig. 1d). This model was solved to steady-state at different spatial positions through the epidermis

Full size image

Images of the epidermis were segmented into discrete tissue layers using histological features apparent with single-target labeling (Fig. 2b; yellow lines) as detailed in our previous work [44]. Next, we introduced a one dimensional ‘layer-normalized’ ordinate which followed the spatiotemporal gradient of keratinocyte terminal differentiation, while minimizing the effect of variations in tissue layer thickness (Fig. 2b) caused by rête ridges (epidermal extensions in to the underlying dermis) and inter-patient differences in the absolute thickness of the epidermis (Additional file 5: Figure S6 in PatALL_-_CompThickness.png). This layer-normalized distance mapped to whole integers at tissue layer boundaries (Table 1), and the relative position within a tissue layer for each sampled region was calculated by linear interpolation, using minimum distances to the deep (d 1 ) and superficial (d 2 ) boundaries of each tissue layer (Fig. 2b). The computational approach to transform fluorescence image data onto a layer-normalized ordinate made the data more comparable across patients (Fig. 2c) and target proteins (Fig. 2d). Fluorescence intensity data reflecting protein abundance were extracted from the sampled pixels and mapped to the layer-normalized distance prior to loess smoothing (Fig. 2c). Smoothed loess curves calculated from these data therefore reflect trends for protein abundance variation across the depth of the epidermis.

Table 1 Spatial partitioning scheme of epidermal tissue layers
Full size table

To use these data for subsequent modeling (Fig. 2d), we mapped normalized signal intensity data from the layer-normalized distance on to discrete spatial partitions which capture the relative number of cells within each tissue layer (Fig. 2c, x-axes; Table 1; further details given in Materials and methods–Image data processing).

Spatial profiles of ERK-MAPK signaling components indicate varying states of the signaling pathway in keratinocytes across the epidermis

For in vitro studies the ERK-MAPK signaling cascade is usually activated by a mitogenic stimulus and changes in protein phosphorylation are examined over minutes-to-hours. Our immunofluorescence data suggest that phosphorylated ERK-MAPK components change in a coordinated manner over the depth of the epidermis (Fig. 3a & b), and thus over a much longer timescale than in vitro studies. The notion that cells within different parts of a tissue behave in a different yet coordinated manner is supported by qualitative observations throughout the literature, thus, we attempted to qualify our observation that ERK-MAPK is modulated across the depth of the epidermis in a quantitative manner.

Fig. 3
figure 3

Spatially coordinated changes to phosphorylated ERK-MAPK components within human epidermis. a Human epidermis simultaneously labeled against phospho-Raf-1 (pS338; cyan), phospho-MEK-1/2 (pS218/pS222; magenta) and phospho-ERK-1/2 (pT183/pY185; yellow). Scale bar represents 10 μm, image data have undergone non-linear transformation to improve printed appearance. b The normalized abundance of cytoplasmic (solid lines) and nuclear (dashed lines) phospho- Raf-1, −MEK-1/2, and ERK-1/2 (colors as above) within interfollicular keratinocytes undergoing terminal differentiation in situ. The c Pearson’s correlation and d mutual information were calculated between all pairwise combinations of target variables using the spatially-conditioned abundance data (grey histograms; axes at left). The strength of these statistical associations was compared to a null distribution calculated from spatially-scrambled data (red probability density function [p.d.f]; axes at right), which was used to calculate two-sided (Pearson’s correlation) and one-sided (mutual information) 99 % confidence intervals (blue vertical lines). For details on the strength of individual relationships please refer to Additional file 6: Table S1

Full size image

Our image processing framework transforms fluorescence intensity data onto the layer-normalized distance ordinate (Fig. 2), producing ‘spatially-conditioned’ data, where the abundance of each target is described with respect to relative spatial position within the epidermis. To test the significance of this spatial conditioning, phospho-protein abundance data were randomly resampled along the spatial axis. These scrambled data were used to establish a null distribution describing the strength of statistical association which could be expected from data with the same marginal distributions (i.e. those arising by chance), in the absence of any spatial information.

We identified the pairwise relationships between phosphorylated ERK-MAPK components which have a spatially-conditioned statistical association which is greater than associations observed within the spatially-scrambled null distribution (Fig. 3c & d; p-value ≤0.01). Many are well-studied, canonical signaling interactions, suggesting that these relationships may be active over the gradient of keratinocyte terminal differentiation (Additional file 6: Table S1). This result highlights the importance of considering spatiotemporal information when examining intracellular signaling mechanisms in situ. Furthermore, given the two week time course associated with keratinocyte terminal differentiation [36–38], this result suggests that phosphorylation of the ERK-MAPK cascade is being regulated over this extended timescale in situ.

Mathematical modeling of ERK-MAPK signaling across the epidermis is consistent with the canonical signaling cascade operating in a quasi-steady-state manner

To investigate whether canonical ERK-MAPK interactions identified in vitro can explain the extended phosphorylation gradient reported here, we use a normalized-Hill differential equation approach [43] to model ERK-MAPK signaling (Fig. 2d). This simplified modeling approach is largely dependent upon the network structure and directionality of interactions defined within the model, making it attractive for our system.

The signaling network (Fig. 1d & Additional file 7: Figure S1 in reaction_network.png) contained the canonical Raf-MEK-ERK interactions with MEK and ERK explicitly separated into cytoplasmic and nuclear compartments with shuttling. Tissue Ca2+ and plasma membrane CaM were included as an activator and inhibitor of cytoplasmic phospho-Raf-1, respectively [18, 23]. Feedback mechanisms within the ERK-MAPK cascade were also incorporated, namely an inhibitory feedback from cytoplasmic phospho-ERK-1/2 to cytoplasmic phospho-Raf-1 [45], and the self dephosphorylation of nuclear phospho-ERK-1/2 through inducing the expression of the nuclear localized ERK phosphatase DUSP4 [46]. State variables are defined as the fractional activation of each signaling species relative to an arbitrary maximal activity (e.g. for ‘ERK’ we examine phospho-ERK-1/2 activity levels).

Relative calcium concentrations (Ca2+) through the depth of the epidermis (Fig. 4a) were derived from the literature [11], and the spatial profile of plasma membrane calmodulin (CaM) (Fig. 4b) was derived from our immunofluorescence data (Additional file 4: Figure S5 in CALM.png). The activity of these input species was used to stimulate the in situ fractional activation of Raf-MEK-ERK using normalized-Hill differential equations. Using spatial steps of ‘half-cell’ width across the epidermis (i.e. steps of 0.5 across 7 spatial partitions; Table 1) we generate conditions for the cell signaling cascade model at 15 spatial locations. For each of these 15 conditions, the model was then integrated until a steady-state was reached. This allows us to determine whether spatial variation in the modeled ‘inputs’ to the signaling cascade can maintain the observed spatial profile of ERK-MAPK activation across the epidermis. Only four parameters were optimized during the least-squares fitting of cytoplasmic and nuclear phospho-ERK-1/2 immunofluorescence data: the baseline levels and amplitudes of Ca2+ and plasma membrane CaM (Table 2).

Fig. 4
figure 4

Simulated and measured abundance profiles of species in our ERK-MAPK signaling model across human epidermis. Our model of the ERK-MAPK pathway is stimulated by: a normalized epidermal Ca2+ as derived by Mauro et al. [11]; and b plasma membrane CaM abundance as derived from our experimental data (Additional file 4: Figure S5i in CALM.png). At each spatial position, the model was run to steady-state with these inputs and c-g the simulated relative abundance of ERK-MAPK components is compared to our in situ experimental measurements. The y-axis in (a) has normalized units, and where experimental data is plotted (b-g), the y-axis represents the z-score of the immunofluorescence pixel intensities. The model input and output abundances were manually scaled to visually fit the experimental data by adjusting the baseline level and amplitude

Full size image
Table 2 Optimized model parameters
Full size table

Qualitatively, the simulated trends for Raf-1, MEK-1/2, and ERK-1/2 fractional activation (Fig. 4c-g) are similar to our spatially-conditioned immunofluorescence data. Within the basal layer (spatial partition 1) the simulations predict relatively low levels of phosphorylation, which is consistent with the immunofluorescence data except for nuclear phospho-MEK-1/2 from two patients (patients 1 and 3). Across the spinous and granular layers (spatial partitions 2–5), there was very good agreement between simulated and measured components, both showing consistent increases in the abundance of the phosphorylated species.

The model predictions performed poorly within the outermost transitional layers (spatial partitions 6 & 7). In particular, the immunofluorescence data showed large decreases in the abundance of cytoplasmic phospho-Raf-1 and both cytoplasmic and nuclear phospho-MEK-1/2, but the model simulated little or no decrease in the phosphorylated abundance of these species. This can probably be attributed to the degradation processes activated during the final stages of keratinocyte differentiation [1], which are not described within our model.

Agreement between our immunofluorescence data and model predictions over the spinous and granular layers suggests that the network structure of the canonical signaling cascade with the specified feedback loops is sufficient to reproduce key features of the extended spatial phosphorylation gradient of ERK-MAPK components observed in human skin. As differentiating keratinocytes progress through the depth of the epidermis and change their expression of signaling receptors and adhesion molecules, the activity of signaling pathways within the keratinocytes shifts, and settles to a new, local steady-state of activation. In our model, we represent tissue microenvironment changes through the relative abundances of tissue Ca2+ and plasma-membrane CaM, and at each position in the epidermis we solve our simple model of the ERK-MAPK cascade within keratinocytes to steady-state. Modulation of ERK-MAPK phosphorylation over these extended timescales may play a central role in the pattern formation processes which establish and maintain epidermal homeostasis.

Discussion

In this study, we developed a novel computational approach to quantify ERK-MAPK activation within human epidermis. Our layer-normalized ordinate incorporates histological landmarks which separate distinct epidermal tissue layers. This allowed us to combine protein abundances measured using single-target labeling into an aggregate data set for mathematical modeling. Previous studies have used the relationship between a keratinocyte’s position within the epidermis and its stage of differentiation, examining protein abundance data (fluorescence signal intensity) with respect to the position within the epidermis, and constructing ‘quantitative spatial profiles’ [47, 48]. Immunohistochemistry has also been applied to measure protein abundances within the skin and to quantify variation across different epidermal layers [49]. This is an important distinction between our work where multiple pattern formation processes maintain complex structures and intracellular signaling states, against studies of in vitro cell signaling where there is only random heterogeneity across the spatial domain.

For cell culture experiments a high degree of intercellular heterogeneity has been noted for ERK signaling, such that the population mean can provide a poor representation of the signal within individual cells [50]. Our statistical analysis (Fig. 3) shows that spatially-dependent changes along the gradient of cellular differentiation account for a large degree of the covariance in the phosphorylation of ERK-MAPK components (i.e. there is a relatively high mutual information between the abundance of successive phosphorylated pathway components). The regulation of signaling component compartmentalization is also thought to be critical for specificity of intracellular signaling pathway responses to different stimuli [26, 51], thus we distinguish the nucleus and cytoplasm when measuring the abundance of active ERK-MAPK components. We believe that methodologies using tissue-specific landmarks with information on subcellular localization, as described here, may be useful to combine data sets and investigate the role of intracellular signaling in the physiological functions of tissues.

Motivated by our statistical analysis, we constructed a model of ERK-MAPK interactions using normalized-Hill differential equations [42, 43]. During model fitting the equations for pathway interactions (e.g. phosphorylation) used default parameters, and only four parameters were optimized: the baseline levels and amplitudes of fractional activation for the input Ca2+ and CaM (Table 2). Given that many of our parameters were poorly constrained, and our signaling network contained several double-phosphorylation steps which can contribute to a non-linear response, the ability of the normalized-Hill approach to qualitatively reproduce in situ immunofluorescence data without extensive parameter-fitting was attractive. A primary limitation of this approach, however, is that it restricts our ability to interpret the modeling results in relation to the mechanistic rate constants of kinetic models. Furthermore, our simplified modeling approach could not examine the intercellular heterogeneity present within the data (Additional file 1: Figure S2, Additional file 2: Figure S3, Additional file 3: Figure S4, Additional file 4: Figure S5 in pRaf1.png; pMEK.png; pERK.png; CALM.png), as discussed in detail below.

The model we developed suggests that the canonical ERK-MAPK signaling cascade operating in a cell-autonomous manner (i.e. without direct cell-cell interactions or long range diffusion of ERK signaling molecules) can reproduce general features of the immunofluorescence data measured in situ, when driven by underlying changes in the surrounding tissue microenvironment. In particular, the relative abundances of tissue Ca2+ and plasma-membrane CaM were implemented as drivers of the model with spatial gradients, suggesting these molecules may be sufficient to establish the extended ERK-MAPK phosphorylation gradient we observed within homeostatic human epidermis. This is largely consistent with observations that increasing extracellular Ca2+ ion concentration promotes keratinocyte terminal differentiation in vitro [13]. It is further supported by the plasma-membrane localization of calmodulin within basal keratinocytes, where it may be down-regulating EGFR and Ras signaling to establish a threshold that prevents proliferation at low doses of growth factors [18]. When CaM-mediated inhibition of Raf was removed from the model, there was an increase in ERK-MAPK pathway activation, and the spatial profile of activation was much more linear, likely reflecting the input Ca2+ gradient (Additional file 8: Figure S7 in SEDML_EpidermalMAPK_varySpatPos_execTimeCourse_noCamRafInhib-output_ScatPlot.png).

Although calcium signaling alone appeared sufficient to drive the observed gradient of ERK-MAPK activation in our model, we acknowledge that this does not contain all possible mechanisms which may influence ERK-MAPK signaling or keratinocyte differentiation. A wide array of signals regulate keratinocyte terminal differentiation [8–12] and mathematical methods have been used to identify dysregulated feedback mechanisms which contribute to a variety of skin pathologies [52–54]. Our model does not attempt to explain the mechanisms that establish and maintain the epidermal Ca2+ gradient [12, 14], in part, because increasingly complex mechanisms are still being identified, such as circadian-rhythm regulated changes in transcript abundance which influence the response of keratinocyte stem cells to Ca2+signaling [55]. Furthermore, due to the ‘single time point’ nature of our in situ image data we also cannot consider intracellular Ca2+ oscillations with epidermal keratinocytes [56, 57].

The model provides a good match to the data across the spinous and granular layers where the gradient of keratinocyte differentiation is more gradual and continuous (Fig. 4c-g). Failure to match the data within the basal layer may be partially attributable to heterogeneity in Ca2+ abundance [12] imposed by asynchronous oscillations, such that small fractions of cells are undergoing different cell processes at any given time-point. Differences in short timescale processes such as this oscillatory behavior may have contributed to variation in the basal layer signaling states that were observed for Patients One and Three (Fig. 4). Rule-based formalisms have given insight towards TGF-β1 regulation of keratinocyte proliferation during wound healing [58] and have been used to test alternative hypotheses for the regulation of basal keratinocyte proliferation [59]. Coupling multiple agent-based models together with tissue-gradient data may provide an approach to further examine the tissue-level effects of intercellular heterogeneity and shorter timescale processes such as the Ca2+ oscillations described above.

An aim of this study was to investigate how the ERK-MAPK cascade extensively studied in vitro may contribute to the maintenance of tissue structure in human epidermis. The pattern of ERK-MAPK signaling along the gradient of keratinocyte terminal differentiation appears to be largely consistent with epidermal biology and the ‘divergent effects’ of ERK-1/2 phosphorylation in regulating cell behavior. For example, transient activation of the ERK-MAPK cascade, typically observed in vitro after a relatively large dose of growth factors or other stimuli, tends to promote mitosis and cellular proliferation [18]. Within the basal layer of the epidermis, we observed relatively low ERK-MAPK activity, while a small number of cells show relatively high phospho-MEK abundance with nuclear localized signal (Additional file 2: Figure S3 in pMEK.png [red arrowheads] & Additional file 9: Figure S8 in FullCellSeg_pMEK_IntDistPDFs_edit.png). Unfortunately, due to the single-target labeling applied for this study, the relative activity of phospho-ERK within these cells cannot be directly examined. Furthermore, due to the relatively low frequency of phospho-MEK bright basal keratinocytes, we have not been able to identify such cells within the simultaneously-labeled tissue samples. Despite this, our observation that the basal epidermis contains a small number of phospho-MEK bright cells is consistent with the notion that basal keratinocytes are undergoing asynchronous division [10], such that at any one time, a small number of cells will show a transient increase in MEK activation just prior to G2/M progression [29]. This small number of proliferative basal keratinocytes with a relatively high nuclear phospho-MEK-1/2 abundance may have contributed to poor agreement between simulated and measured abundances (Fig. 4e).

As keratinocytes detach from the underlying basement membrane and begin execution of their terminal differentiation program, there is a small stepwise increase in the abundance of phospho-MEK-1/2 followed by a sustained, gradual increase over the spinous and granular tissue layers. Increased phosphorylation of ERK-1/2 has been observed with terminal differentiation in various cell lines [18, 60] and sustained ERK-MAPK pathway activation has been reported to induce growth arrest in epidermoid carcinoma cells [19]. It is tempting to speculate that sustained activation of the ERK-MAPK pathway helps to suppress cellular proliferation and ensure growth arrest in suprabasal keratinocytes by reducing the availability of unphosphorylated MEK-1/2 and ERK-1/2, preventing the transient increases that promote mitosis. The presence of complex regulatory processes described above, and the abnormal differentiation programs adopted by keratinocytes in vitro [39–41], make it difficult to test model predictions directly; however, transgenic mice with modified signaling components under inducible expression within the epidermis may be useful for testing such mechanisms [17, 61]. We believe that the observations in this work highlight the need for further studies of in situ signaling pathway activation.

Conclusions

This study developed a novel, joint experimental-computational framework for investigating the regulation of ERK-MAPK signaling in homeostatic human epidermis. We have identified an extended phospho-protein gradient of ERK-MAPK signaling cascade components over the depth of the epidermis, and we utilized statistical analysis and mathematical modeling to demonstrate that canonical ERK-MAPK interactions driven by Ca2+ signaling components operating at quasi-steady-state are sufficient to establish this extended phosphorylation gradient. These findings may have implications for the wide range of cancer drugs which therapeutically target ERK-MAPK signaling components.

Materials and methods

Ethics statement

Human skin samples were obtained with written informed consent under a protocol approved by the New Zealand Northern Regional X Ethics Committee (project number NTX/08/09/086) and Counties-Manukau District Health Board (project number 681).

Tissue collection

Fresh skin samples obtained from healthy patients undergoing plastic or reconstructive surgery were snap frozen in liquid nitrogen and stored at −80 °C.

Immunofluorescence labeling and confocal microscopy

A detailed description of the immunofluorescence labeling and confocal microscopy protocols, including antibody catalogue numbers, are given elsewhere [44]. The antibody against calmodulin (UniProt: P62158) targets multiple isoforms. A phospho-serine 338 specific antibody was used to examine the activation state of Raf-1 (UniProt: P04049), as this phospho-epitope is routinely used as a surrogate marker for Raf-1 activation [45]. Antibodies against double-phosphorylated MEK-1/2 (pS218/pS222) and double-phosphorylated ERK-1/2 (pT183/pY185) were used to examine the abundance of active MEK-1/2 (UniProt: Q02750/P36507) and ERK-1/2 (UniProt: P27361/P28482).

Human epidermis was sectioned within a cryostat at −20 °C and immediately fixed with paraformaldehyde at room temperature, then labeled by indirect immunofluorescence and imaged using a Leica TCS SP2 confocal microscope. Secondary antibody control slides were imaged at the same gain and offset settings to ensure that the effects of non-specific labeling or auto-fluorescence were minimal. An imaging resolution at or exceeding the Nyquist sampling criteria was used, allowing the plasma-membrane, nucleus and cytoplasm of individual cells to be identified. Gain and offset values were adjusted to minimize underflow and overflow, and under these imaging conditions, fluorescence signal intensity was interpreted as an approximately linear measure of protein abundance within each voxel [62]. Image data were stored in an uncompressed .TIFF format.

Image data processing

A collection of scripts for quantifying the image data was written using MATLAB (MathWorks, Natick MA). A graphical user interface was developed using the MATLAB image processing toolbox and the ginput2 function (MATLAB Central File Exchange FileID: #20645) to facilitate selection of regions of interest (ROIs). Tissue layer boundaries were specified using the path tool within GIMP (v. 2.4.7) and exported as a binarized .TIFF format. The minimum distance between ROIs and surrounding tissue layer boundaries was determined using the intrinsic MATLAB function dist. For analyses where data were combined across patients (Fig. 3), the signal intensity data were z-score normalized for each patient and then combined.

Loess smoothing was chosen due to its ability to cope with edge effects (by selecting local sample points dependent upon the smoothing frequency parameter) [63]. Loess smoothing was applied to the signal intensity data across each tissue layer using the intrinsic MATLAB function smooth. Although the normalized distance co-ordinate system captured the distinct epidermal layers well, signal intensity data examined against this gradient did not capture the different thicknesses of each tissue layer (Additional file 5: Figure S6 in PatALL_-_CompThickness.png). Thus, the spatial domain was further partitioned for discrete spatial sampling, in a manner that reflected the ratio of 1:4:2 cells between the tissue layers (Table 1).

Statistical test of spatial conditioning

Image data (Fig. 3; Additional file 1: Figure S2, Additional file 2: Figure S3, Additional file 3: Figure S4) suggested that activation of the ERK-MAPK cascade occurs in a coordinated manner over the spatiotemporal gradient of keratinocyte terminal differentiation. To test this in a quantitative manner, a null hypothesis was proposed:

“the strength of statistical association between measured fluorescence signals (signaling component abundance) is independent of spatial conditioning”

Signal intensity data were extracted over 28 spatial partitions (Table 1). Data were examined over the spatial dimension, and the Pearson’s correlation and mutual information between each signaling component was calculated using the MATLAB corr and mutualinfo (MATLAB Central File Exchange FileID: #14888) functions, respectively. Data scrambling was applied to abrogate spatial conditioning of each variable while retaining the same marginal distributions of signal intensity data. For each target, signal intensity data were scrambled along the spatial domain (i.e. random sampling without replacement) 100,000 times. From these scrambled data, the Pearson’s correlation and mutual information were calculated between targets and the probability density functions were taken as null distributions. The two-sided (Pearson’s correlation) and one-sided (mutual information) 99 % confidence intervals were identified as thresholds (Fig. 3c & d; blue vertical lines). For the spatially-conditioned data, relationships with a statistical association exceeding these thresholds (Additional file 6: Table S1) were considered as evidence against the null hypothesis that the statistical associations are independent of spatial conditioning (p ≤ 0.01).

Estimating cytoplasmic-to-nuclear volume ratio

Scripts were written using MATLAB to manually input sample points that demarcate the edge of cellular domains at various z-positions through the image data. A series of image processing functions were applied to: connect these data points and fill the selected area, smooth edges and subtract the nucleus from the cytoplasm, and interpolate the segmented areas between z-positions to produce a volumetric segmentation. From these data, the volume ratio of cytoplasm-to-nucleus was calculated for a number of cells across the depth of the epidermis, as shown in Additional file 10: Figure S10 in cyto_nuc_ratios_edit.png.

Mathematical modeling

An SBML [64] implementation of the model has been deposited in BioModels Database [65] and assigned the identifier MODEL1503270000. Scripts used in this project, including SED-ML [66] to execute the model at different spatial positions and reproduce results from this paper are available from http://sourceforge.net/projects/EpidermalERKMAPK/.

The reaction kinetic diagram describing the signaling processes (Additional file 7: Figure S1 in reaction_network.png), and the simplified version of the network used for the model (Fig. 1d) were drawn as SBGN [67] process diagrams using VANTED [68]. To prevent over-parameterization of the model, indirect signaling processes that involve unmeasured components were ‘collapsed’ into a single arc. For example: Ca2+ is a direct input to phospho-Raf-1 within the model, as intermediate variables such as Ras-GTP were not measured; and similarly the role of nuclear phospho-ERK-1/2 in promoting its own dephosphorylation is modeled as a self-interaction, as the induced phosphatase DUSP4 was not measured.

Three types of interactions between species were used to construct the mathematical model:

  1. 1

    Activation: Represented with a normalized activating Hill function [43]:

    $$ {f}_{act}(X)=\frac{B{X}^n}{K^n+{X}^n} $$
    (1.1)

    where B and K are constrained such that fact(0) = 0, fact(EC50) = 0.5 and fact(1) = 1. EC50 is the fractional activation of an input species at which a half-maximal activation of an output species is induced;

  2. 2.

    Inhibition: Represented as the negative of activation;

  3. 3.

    Dephosphorylation: Represented as the product of inhibition with the fractional activation of the output species, as the abundance of the output species determines the maximal level of dephosphorylation that can occur.

In addition, shuttling between the cytoplasmic and nuclear compartments was modeled to be linearly proportional to the fractional activations of the species in the respective compartments. Pathway crosstalk was implemented by summing individual species interactions. Using this approach, the ERK-MAPK signaling network in Fig. 1d was represented with the following differential equations:

$$ \begin{array}{c}\hfill \frac{dRafc}{dt}=\frac{1}{\tau_{Rafc}}\left( \max \left({w}_{Ca, Rafc}{f}_{act}(Ca)-{w}_{CaM, Rafc}{f}_{act}(CaMm)-{w}_{ERKc, Rafc}{f}_{act}(ERKc),0\right)Raf{c}_{\max }- Rafc\right)\hfill \\ {}\hfill \frac{dMEKc}{dt}=\frac{1}{\tau_{MEKc}}\left({w}_{Rafc, MEKc}{f}_{act}(Rafc)MEK{c}_{\max }- MEKc-{s}_{MEKc n} MEKc+\frac{s_{MEKn c}}{r} MEKn\right)\hfill \\ {}\hfill \frac{dMEKn}{dt}=\frac{1}{\tau_{MEKn}}\left(- MEKn-{s}_{MEKn c} MEKn+r{s}_{MEKc n} MEKc\right)\hfill \\ {}\hfill \frac{dERKc}{dt}=\frac{1}{\tau_{ERKc}}\left({w}_{MEKc, ERKc}{f}_{act}(MEKc)ERK{c}_{\max }- ERKc-{s}_{ERKc n} ERKc+\frac{s_{ERKn c}}{r} ERKn\right)\hfill \\ {}\hfill \frac{dERKn}{dt}=\begin{array}{c}\hfill \frac{1}{\tau_{ERKn}}\Big(\left({w}_{MEKn, ERKn}{f}_{act}(MEKn)-{w}_{ERKn, ERKn}{f}_{act}(ERKn) ERKn\right)ERK{n}_{\max }- ERKn\hfill \\ {}\hfill -{s}_{ERKn c} ERKn+r{s}_{ERKc n} ERKc\Big)\hfill \end{array}\hfill \end{array} $$
(1.2)

where: state variables Ca, CaM, Rafc, MEKc, MEKn, ERKc and ERKn represent fractional activation of tissue Ca2+, plasma membrane CaM, cytoplasmic phospho-Raf-1, cytoplasmic phospho-MEK-1/2, nuclear phospho-MEK-1/2, cytoplasmic phospho-ERK-1/2, and nuclear phospho-ERK-1/2, respectively (note only active/phosphorylated ERK-MAPK components were examined in this study–although immunofluorescence data were collected for the non-phosphorylated forms, technical issues for MEK in particular prevented their inclusion within the quantitative modeling, as illustrated in Additional file 11: Figure S9 in MEK.png); τ is the time constant for a given species; w is the reaction weight; X max is the maximal fractional activation of species; s Xcn and s Xnc are the cytoplasmic-nuclear and nuclear-cytoplasmic shuttling parameters for species X, respectively; and r is the cytoplasmic-nuclear volume ratio. In the equation governing cytoplasmic phospho-Raf-1 activity, the max function was incorporated to enforce that the inhibitory effects of plasma membrane CaM and cytoplasmic phospho-ERK-1/2 do not escalate into dephosphorylation effects instead (i.e. these inhibitions can only negate the activating effect of extracellular Ca2+, but cannot reduce fractional activation of cytoplasmic phospho-Raf-1). The scaling by r of the shuttling terms was to account for volume effects between the cytoplasm and nucleus. The cytoplasmic-nuclear volume ratio was set to vary linearly from 2 (deep) to 5 (superficial), as derived from the immunofluorescence data (Additional file 10: Figure S10 in cyto_nuc_ratios_edit.png).

Default parameters were used for the normalized Hill differential equations (n = 1.4, EC50 = 0.5, τ = 1, w = 1 and Xmax = 1) [43]. MEK and ERK shuttling parameters (sMEKcn = 0.05, sMEKnc = 0.5, sERKcn = 0.01 and sERKnc = 0.01) were selected to reflect experimentally-determined ratios of cytoplasm-to-nucleus and nucleus-to-cytoplasm transport coefficients [27]. It is difficult to directly reconcile the parameter values derived for reaction kinetic models of individual cells undergoing rapid transient changes, and those used here for examining quasi-steady state changes over an extended timescale. One would expect similar shuttling rates associated with biochemical features, such as the nuclear export sequence on MEK [69]. There may be differences, however, as the HeLa cells used by Fujioka et al. [27] are derived from an adenocarcinoma and thus likely express different levels of scaffolding molecules/co-factors which vary between tissues [25, 70]. These issues highlight the need for studies which examine these signaling pathways within normal cells and homeostatic tissues.

The cell-autonomous signaling model (above) is driven by spatial profiles of extracellular Ca2+ and plasma membrane CaM activation levels across the epidermis. The Ca2+profile was modeled based on a previous experimental report [11] and CaM based on our own immunofluorescence data respectively as:

$$ \begin{array}{c}\hfill Ca={b}_{Ca}+{a}_{Ca}\left\{\begin{array}{ccc}\hfill d/5\hfill & \hfill for\hfill & \hfill d\in \left[0,5\right)\hfill \\ {}\hfill \left(7-d\right)/2\hfill & \hfill for\hfill & \hfill d\in \left[5,7\right]\hfill \end{array}\right.\hfill \\ {}\hfill CaMm={b}_{CaM}+{a}_{CaM}\left\{\begin{array}{ccc}\hfill \left(d+1\right)/2\hfill & \hfill for\hfill & \hfill d\in \left[0,1\right)\hfill \\ {}\hfill {e}^{\left(1-d\right)}\hfill & \hfill for\hfill & \hfill d\in \left[1,5\right)\hfill \\ {}\hfill 0\hfill & \hfill for\hfill & \hfill d\in \left[5,7\right]\hfill \end{array}\right.\hfill \end{array} $$
(1.3)

where: b and a are the baseline levels and amplitudes of the input fractional activations respectively (Table 2); d specifies the spatial partition (i.e. depth) within the epidermis.

A least-squares fit to the aggregate cytoplasmic and nuclear phospho-ERK-1/2 immunofluorescence data was performed with the MATLAB lsqnonlin and polyfit functions to optimize values of b and a. Differential equations were implemented in MATLAB (v7.12; R2011b) and solved numerically for each spatial partition using the ode23 function to find the steady-state (50 time units, with a time step of 0.1 time units).

References

  1. Candi E, Schmidt R, Melino G. The cornified envelope: a model of cell death in the skin. Nat Rev Mol Cell Biol. 2005;6:328–40.

    Article  CAS  PubMed  Google Scholar 

  2. Lippens S, Denecker G, Ovaere P, Vandenabeele P, Declercq W. Death penalty for keratinocytes: apoptosis versus cornification. Cell Death Differ. 2005;12 Suppl 2:1497–508.

    Article  CAS  PubMed  Google Scholar 

  3. Adams JC, Watt FM. Fibronectin inhibits the terminal differentiation of human keratinocytes. Nature. 1989;340:307–9.

    Article  CAS  PubMed  Google Scholar 

  4. Fuchs E. Epidermal differentiation: the bare essentials. J Cell Biol. 1990;111:2807–14.

    Article  CAS  PubMed  Google Scholar 

  5. Eckert RL, Crish JF, Robinson NA. The epidermal keratinocyte as a model for the study of gene regulation and cell differentiation. Physiol Rev. 1997;77:397–424.

    CAS  PubMed  Google Scholar 

  6. Angel P, Szabowski A, Schorpp-Kistner M. Function and regulation of AP-1 subunits in skin physiology and pathology. Oncogene. 2001;20:2413–23.

    Article  CAS  PubMed  Google Scholar 

  7. Blanpain C, Fuchs E. Epidermal homeostasis: a balancing act of stem cells in the skin. Nat Rev Mol Cell Biol. 2009;10:207–17.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  8. Fuchs E, Dowling J, Segre J, Lo SH, Yu QC. Integrators of epidermal growth and differentiation: distinct functions for beta 1 and beta 4 integrins. Curr Opin Genet Dev. 1997;7:672–82.

    Article  CAS  PubMed  Google Scholar 

  9. Sonnenberg A, Calafat J, Janssen H, Daams H, van der Raaij-Helmer LM, Falcioni R, et al. Integrin alpha 6/beta 4 complex is located in hemidesmosomes, suggesting a major role in epidermal cell-basement membrane adhesion. J Cell Biol. 1991;113:907–17.

    Article  CAS  PubMed  Google Scholar 

  10. Muller EJ, Williamson L, Kolly C, Suter MM. Outside-in signaling through integrins and cadherins: a central mechanism to control epidermal growth and differentiation? J Invest Dermatol. 2008;128:501–16.

    Article  CAS  PubMed  Google Scholar 

  11. Mauro T, Bench G, Sidderas-Haddad E, Feingold K, Elias P, Cullander C. Acute barrier perturbation abolishes the Ca2+ and K+ gradients in murine epidermis: quantitative measurement using PIXE. J Invest Dermatol. 1998;111:1198–201.

    Article  CAS  PubMed  Google Scholar 

  12. Celli A, Sanchez S, Behne M, Hazlett T, Gratton E, Mauro T. The epidermal Ca(2+) gradient: measurement using the phasor representation of fluorescent lifetime imaging. celli.anna@gmail.com. Biophys J. 2010;98:911–21.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  13. Yuspa SH, Kilkenny AE, Steinert PM, Roop DR. Expression of murine epidermal differentiation markers is tightly regulated by restricted extracellular calcium concentrations in vitro. J Cell Biol. 1989;109:1207–17.

    Article  CAS  PubMed  Google Scholar 

  14. Elias P, Ahn S, Brown B, Crumrine D, Feingold KR. Origin of the epidermal calcium gradient: regulation by barrier status and role of active vs passive mechanisms. J Invest Dermatol. 2002;119:1269–74.

    Article  CAS  PubMed  Google Scholar 

  15. Iversen L, Johansen C, Kragballe K. Signal transduction pathways in human epidermis. Eur J Dermatol. 2005;15:4–12.

    CAS  PubMed  Google Scholar 

  16. Eckert RL, Efimova T, Dashti SR, Balasubramanian S, Deucher A, Crish JF, et al. Keratinocyte survival, differentiation, and death: many roads lead to mitogen-activated protein kinase. J Investig Dermatol Symp Proc. 2002;7:36–40.

    Article  CAS  PubMed  Google Scholar 

  17. Scholl FA, Dumesic PA, Khavari PA. Mek1 alters epidermal growth and differentiation. Cancer Res. 2004;64:6035–40.

    Article  CAS  PubMed  Google Scholar 

  18. Agell N, Bachs O, Rocamora N, Villalonga P. Modulation of the Ras/Raf/MEK/ERK pathway by Ca(2+), and calmodulin. Cell Signal. 2002;14:649–54.

    Article  CAS  PubMed  Google Scholar 

  19. Bromberg JF, Fan Z, Brown C, Mendelsohn J, Darnell JE. Epidermal growth factor-induced growth inhibition requires Stat1 activation. Cell Growth Differ. 1998;9:505–12.

    CAS  PubMed  Google Scholar 

  20. Jost M, Huggett TM, Kari C, Rodeck U. Matrix-independent survival of human keratinocytes through an EGF receptor/MAPK-kinase-dependent pathway. Mol Biol Cell. 2001;12:1519–27.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  21. Manohar A, Shome SG, Lamar J, Stirling L, Iyer V, Pumiglia K, et al. Alpha 3 beta 1 integrin promotes keratinocyte cell survival through activation of a MEK/ERK signaling pathway. J Cell Sci. 2004;117:4043–54.

    Article  CAS  PubMed  Google Scholar 

  22. Mainiero F, Murgia C, Wary KK, Curatola AM, Pepe A, Blumemberg M, et al. The coupling of α6β4 integrin to Ras–MAP kinase pathways mediated by Shc controls keratinocyte proliferation. EMBO J. 1997;16:2365–75.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  23. Schmidt M, Goebeler M, Posern G, Feller SM, Seitz CS, Brocker EB, et al. Ras-independent activation of the Raf/MEK/ERK pathway upon calcium-induced differentiation of keratinocytes. J Biol Chem. 2000;275:41011–7.

    Article  CAS  PubMed  Google Scholar 

  24. Eferl R, Wagner EF. AP-1: a double-edged sword in tumorigenesis. Nat Rev Cancer. 2003;3:859–68.

    Article  CAS  PubMed  Google Scholar 

  25. Kolch W. Coordinating ERK/MAPK signalling through scaffolds and inhibitors. Nat Rev Mol Cell Biol. 2005;6:827–37.

    Article  CAS  PubMed  Google Scholar 

  26. Pouyssegur J, Volmat V, Lenormand P. Fidelity and spatio-temporal control in MAP kinase (ERKs) signalling. Biochem Pharmacol. 2002;64:755–63.

    Article  CAS  PubMed  Google Scholar 

  27. Fujioka A, Terai K, Itoh RE, Aoki K, Nakamura T, Kuroda S, et al. Dynamics of the Ras/ERK MAPK cascade as monitored by fluorescent probes. J Biol Chem. 2006;281:8917–26.

    Article  CAS  PubMed  Google Scholar 

  28. Huang W, Alessandrini A, Crews CM, Erikson RL. Raf-1 forms a stable complex with Mek1 and activates Mek1 by serine phosphorylation. Proc Natl Acad Sci U S A. 1993;90:10947–51.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  29. Harding A, Giles N, Burgess A, Hancock JF, Gabrielli BG. Mechanism of mitosis-specific activation of MEK1. J Biol Chem. 2003;278:16747–54.

    Article  CAS  PubMed  Google Scholar 

  30. Bitangcol JC, Chau AS, Stadnick E, Lohka MJ, Dicken B, Shibuya EK. Activation of the p42 mitogen-activated protein kinase pathway inhibits Cdc2 activation and entry into M-phase in cycling Xenopus egg extracts. Mol Biol Cell. 1998;9:451–67.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  31. Kolch W, Calder M, Gilbert D. When kinases meet mathematics: the systems biology of MAPK signalling. FEBS Lett. 2005;579:1891–5.

    Article  CAS  PubMed  Google Scholar 

  32. Schoeberl B, Eichler-Jonsson C, Gilles ED, Müller G. Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors. Nat Biotechnol. 2002;20:370–5.

    Article  PubMed  Google Scholar 

  33. Hornberg JJ, Binder B, Bruggeman FJ, Schoeberl B, Heinrich R, Westerhoff HV. Control of MAPK signalling: from complexity to what really matters. Oncogene. 2005;24:5533–42.

    Article  CAS  PubMed  Google Scholar 

  34. Reeves GT, Muratov CB, Schupbach T, Shvartsman SY. Quantitative models of developmental pattern formation. Dev Cell. 2006;11:289–300.

    Article  CAS  PubMed  Google Scholar 

  35. Kim Y, Coppey M, Grossman R, Ajuria L, Jimenez G, Paroush Z, et al. MAPK substrate competition integrates patterning signals in the Drosophila embryo. Curr Biol. 2010;20:446–51.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  36. Weinstein GD, Scott EJV. Autoradiographic analysis of turnover times of normal and psoriatic epidermis. J Invest Dermatol. 1965;45:257–62.

    Article  CAS  PubMed  Google Scholar 

  37. Halprin KM. Epidermal “turnover time”–a re-examination. Br J Dermatol. 1972;86:14–9.

    Article  CAS  PubMed  Google Scholar 

  38. Potten CS, Saffhill R, Maibach HI. Measurement of the transit time for cells through the epidermis and stratum corneum of the mouse and guinea-pig. Cell Tissue Kinet. 1987;20:461–72.

    CAS  PubMed  Google Scholar 

  39. Eichner R, Bonitz P, Sun TT. Classification of epidermal keratins according to their immunoreactivity, isoelectric point, and mode of expression. J Cell Biol. 1984;98:1388–96.

    Article  CAS  PubMed  Google Scholar 

  40. Weiss RA, Eichner R, Sun TT. Monoclonal antibody analysis of keratin expression in epidermal diseases: a 48- and 56-kdalton keratin as molecular markers for hyperproliferative keratinocytes. J Cell Biol. 1984;98:1397–406.

    Article  CAS  PubMed  Google Scholar 

  41. Olsen E, Rasmussen HH, Celis JE. Identification of proteins that are abnormally regulated in differentiated cultured human keratinocytes. Electrophoresis. 1995;16:2241–8.

    Article  CAS  PubMed  Google Scholar 

  42. Ryall KA, Holland DO, Delaney KA, Kraeutler MJ, Parker AJ, Saucerman JJ. Network reconstruction and systems analysis of cardiac myocyte hypertrophy signaling. J Biol Chem. 2012;287:42259–68.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  43. Kraeutler MJ, Soltis AR, Saucerman JJ. Modeling cardiac β-adrenergic signaling with normalized-Hill differential equations: comparison with a biochemical model. BMC Syst Biol. 2010;4:157.

    Article  PubMed Central  PubMed  Google Scholar 

  44. Cursons J, Hurley D, Angel CE, Dunbar R, Crampin EJ, Jacobs MD. Inference of an in situ epidermal intracellular signaling cascade. Conf Proc IEEE Eng Med Biol Soc. 2010;2010:799–802.

    PubMed  Google Scholar 

  45. Matallanas D, Birtwistle M, Romano D, Zebisch A, Rauch J, von Kriegsheim A, et al. Raf family kinases: old dogs have learned new tricks. Genes Cancer. 2011;2:232–60.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  46. Cagnol S, Rivard N. Oncogenic KRAS and BRAF activation of the MEK/ERK signaling pathway promotes expression of dual-specificity phosphatase 4 (DUSP4/MKP2) resulting in nuclear ERK1/2 inhibition. Oncogene. 2012;32:564–76.

    Article  PubMed  Google Scholar 

  47. Grabe N, Pommerencke T, Steinberg T, Dickhaus H, Tomakidi P. Reconstructing protein networks of epithelial differentiation from histological sections. Bioinformatics. 2007;23:3200–8.

    Article  CAS  PubMed  Google Scholar 

  48. Pommerencke T, Steinberg T, Dickhaus H, Tomakidi P, Grabe N. Nuclear staining and relative distance for quantifying epidermal differentiation in biomarker expression profiling. BMC Bioinformatics. 2008;9:473.

    Article  PubMed Central  PubMed  Google Scholar 

  49. Kokolakis G, Panagis L, Stathopoulos E, Giannikaki E, Tosca A, Krüger-Krasagakis S. From the protein to the graph: how to quantify immunohistochemistry staining of the skin using digital imaging. J Immunol Methods. 2008;331:140–6.

    Article  CAS  PubMed  Google Scholar 

  50. Walker DC, Georgopoulos NT, Southgate J. From pathway to population–a multiscale model of juxtacrine EGFR-MAPK signalling. BMC Syst Biol. 2008;2:102.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  51. Lipshtat A, Neves SR, Iyengar R. Specification of spatial relationships in directed graphs of cell signaling networks. Ann N Y Acad Sci. 2009;1158:44–56.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  52. Tanaka RJ, Ono M, Harrington HA. Skin barrier homeostasis in atopic dermatitis: feedback regulation of kallikrein activity. PLoS One. 2011;6:e19895.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  53. Tanaka RJ, Ono M. Skin disease modeling from a mathematical perspective. J Invest Dermatol. 2013;133:1472–8.

    Article  CAS  PubMed  Google Scholar 

  54. Dominguez-Huttinger E, Ono M, Barahona M, Tanaka RJ. Risk factor-dependent dynamics of atopic dermatitis: modelling multi-scale regulation of epithelium homeostasis. Interface Focus. 2013;3:20120090.

    Article  PubMed Central  PubMed  Google Scholar 

  55. Janich P, Toufighi K, Solanas G, Luis NM, Minkwitz S, Serrano L, et al. Human epidermal stem cell function is regulated by circadian oscillations. Cell Stem Cell. 2013;13:745–53.

    Article  CAS  PubMed  Google Scholar 

  56. Denda M, Denda S. Air-exposed keratinocytes exhibited intracellular calcium oscillation. Skin Res Technol. 2007;13:195–201.

    Article  PubMed  Google Scholar 

  57. Tsutsumi M, Inoue K, Denda S, Ikeyama K, Goto M, Denda M. Mechanical-stimulation-evoked calcium waves in proliferating and differentiated human keratinocytes. Cell Tissue Res. 2009;338:99–106.

    Article  CAS  PubMed  Google Scholar 

  58. Adra S, Sun T, MacNeil S, Holcombe M, Smallwood R. Development of a three dimensional multiscale computational model of the human epidermis. PLoS One. 2010;5:e8511.

    Article  PubMed Central  PubMed  Google Scholar 

  59. Li X, Upadhyay AK, Bullock AJ, Dicolandrea T, Xu J, Binder RL, et al. Skin stem cell hypotheses and long term clone survival–explored using agent-based modelling. Sci Rep. 2013;3:1904.

    CAS  PubMed Central  PubMed  Google Scholar 

  60. Adachi T, Kar S, Wang M, Carr BI. Transient and sustained ERK phosphorylation and nuclear translocation in growth control. J Cell Physiol. 2002;192:151–9.

    Article  CAS  PubMed  Google Scholar 

  61. Tarutani M, Cai T, Dajee M, Khavari PA. Inducible activation of Ras and Raf in adult epidermis. Cancer Res. 2003;63:319–23.

    CAS  PubMed  Google Scholar 

  62. Murray JM. Practical aspects of quantitative confocal microscopy. Methods Cell Biol. 2007;81:467–78.

    Article  PubMed  Google Scholar 

  63. Cleveland WS, Devlin SJ. Locally weighted regression: an approach to regression analysis by local fitting. J Am Stat Assoc. 1988;83:596–610.

    Article  Google Scholar 

  64. Hucka M, Finney A, Sauro HM, Bolouri H, Doyle JC, Kitano H, et al. The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics. 2003;19:524–31.

    Article  CAS  PubMed  Google Scholar 

  65. Li C, Donizelli M, Rodriguez N, Dharuri H, Endler L, Chelliah V, et al. BioModels database: an enhanced, curated and annotated resource for published quantitative kinetic models. BMC Syst Biol. 2010;4:92.

    Article  PubMed Central  PubMed  Google Scholar 

  66. Waltemath D, Adams R, Bergmann FT, Hucka M, Kolpakov F, Miller AK, et al. Reproducible computational biology experiments with SED-ML–the Simulation Experiment Description Markup Language. BMC Syst Biol. 2011;5:198.

    Article  PubMed Central  PubMed  Google Scholar 

  67. Le Novere N, Hucka M, Mi H, Moodie S, Schreiber F, Sorokin A, et al. The systems biology graphical notation. Nat Biotechnol. 2009;27:735–41.

    Article  PubMed  Google Scholar 

  68. Junker BH, Klukas C, Schreiber F. VANTED: a system for advanced data analysis and visualization in the context of biological networks. BMC Bioinformatics. 2006;7:109.

    Article  PubMed Central  PubMed  Google Scholar 

  69. Fukuda M, Gotoh I, Adachi M, Gotoh Y, Nishida E. A novel regulatory mechanism in the mitogen-activated protein (MAP) kinase cascade. Role of nuclear export signal of MAP kinase kinase. J Biol Chem. 1997;272:32642–8.

    Article  CAS  PubMed  Google Scholar 

  70. Witzel F, Maddison L, Bluthgen N. How scaffolds shape MAPK signaling: what we know and opportunities for systems approaches. Front Physiol. 2012;3:475.

    Article  PubMed Central  PubMed  Google Scholar 

Download references

Acknowledgements

We gratefully acknowledge the contributions of Jenni Chen and Dr Kate Angel during development of the immunofluorescence labeling methods. We would also like to thank Dr Melissa Davis for critical proof-reading of the manuscript.

JC received funding from the Tertiary Education Commission (NZ) Top Achiever Doctoral Scholarship. This project received funding through a New Economy Research Fund (NZ) grant for Advanced Skin Imaging (contract C08X0801). This research was in part conducted and funded by the Australian Research Council Centre of Excellence in Convergent Bio-Nano Science and Technology (project number CE140100036).

Author information

Authors and Affiliations

  1. Systems Biology Laboratory, Melbourne School of Engineering, University of Melbourne, Melbourne, Australia

    Joseph Cursons, Jerry Gao, Daniel G. Hurley & Edmund J. Crampin

  2. NICTA Victoria Research Lab, Melbourne, Australia

    Joseph Cursons & Daniel G. Hurley

  3. ARC Centre of Excellence in Convergent Bio-Nano Science and Technology, Melbourne School of Engineering, University of Melbourne, Melbourne, Australia

    Joseph Cursons & Edmund J. Crampin

  4. Auckland Bioengineering Institute, University of Auckland, Auckland, New Zealand

    Joseph Cursons, Daniel G. Hurley & Edmund J. Crampin

  5. Maurice Wilkins Centre, University of Auckland, Auckland, New Zealand

    Joseph Cursons, Daniel G. Hurley, Cristin G. Print, P. Rod Dunbar & Edmund J. Crampin

  6. Bioinformatics Institute, University of Auckland, Auckland, New Zealand

    Daniel G. Hurley & Cristin G. Print

  7. Faculty of Medical and Health Sciences, University of Auckland, Auckland, New Zealand

    Daniel G. Hurley & Cristin G. Print

  8. School of Biological Sciences, University of Auckland, Auckland, New Zealand

    P. Rod Dunbar

  9. Department of Biology, New Zealand International College, ACG New Zealand, Auckland, New Zealand

    Marc D. Jacobs

  10. School of Mathematics and Statistics, University of Melbourne, Melbourne, Australia

    Edmund J. Crampin

  11. School of Medicine, University of Melbourne, Melbourne, Australia

    Edmund J. Crampin

Authors
  1. Joseph Cursons
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Jerry Gao
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Daniel G. Hurley
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Cristin G. Print
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. P. Rod Dunbar
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Marc D. Jacobs
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Edmund J. Crampin
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Edmund J. Crampin.

Additional information

Competing interests

The author's declare that they have no competing interests.

Authors’ contributions

JC planned the experiments, performed immunofluorescence labeling and confocal microscopy imaging, developed computational tools for processing and analyzing the image data, produced the SBML representation of the model and SED-ML scripts to reproduce results shown in the manuscript, and wrote the manuscript. JG developed the mathematical model and MATLAB code, and wrote the manuscript. DH helped develop computational tools for data analysis and proof-read the manuscript. CP, PRD and MDJ planned the experiments and assisted with data analysis. EJC planned the experiments, assisted with data analysis, helped design the mathematical model and wrote the manuscript. All authors read and approved the final manuscript.

Additional files

Additional file 1: Figure S2.

Human epidermis (Patient Two) labeled against phospho-Raf-1 (pS338) using immunofluorescence labeling with confocal microscopy imaging. (a, b) Immunofluorescence images are displayed together with (c) a surface rendering of the signal intensity within a suprabasal keratinocyte. The z-score normalized sampled signal intensity data (data points; plotted relative to the sample mean, \( \overline{\mu} \); and sample standard deviation, \( \overline{\sigma} \)) and loess smoothed signals (solid lines) associated with the (d) cytoplasm and (f) nuclei are displayed for Patient One (red), Two (green) and Three (blue), together with the 90 % confidence interval for positive and negative residuals (dashed lines), and the sampled data clouds for Patient Two. Histograms of the discretized signal intensity within the (e) cytoplasm and (g) nucleus across individual spatial divisions are shown, together with the associated loess-curves (magenta line). The regions displayed in (b) and (c) are highlighted within (a) by the white and orange dashed lines, respectively. Image data have undergone a non-linear transformation to improve visual appearance.

Additional file 2: Figure S3.

Human epidermis (Patient Two) labeled against phospho-MEK-1/2 (pS218/pS222) using immunofluorescence labeling with confocal microscopy imaging. (a, b) Immunofluorescence images are displayed together with (c) a surface rendering of the signal intensity within suprabasal keratinocytes). The z-score normalized sampled signal intensity data (data points; plotted relative to the sample mean, \( \overline{\mu} \); and sample standard deviation, \( \overline{\sigma} \)) and loess smoothed signals (solid lines) associated with the (d) cytoplasm and (f) nuclei are displayed for Patient One (red), Two (green) and Three (blue), together with the 90 % confidence interval for positive and negative residuals (dashed lines), and the sampled data clouds for Patient Two. Histograms of the discretized signal intensity within the (e) cytoplasm and (g) nuclei across individual spatial divisions are shown, together with the associated loess-curves (magenta line). The regions displayed in (b) and (c) are highlighted within (a) by the white and orange dashed lines, respectively. Scale bars represent 10 μm. Image data have undergone a non-linear transformation to improve visual appearance. Basal cells with a relatively-high phospho-MEK-1/2 signal intensity are highlighted (red arrowhead).

Additional file 3: Figure S4.

Human epidermis (Patient One) labeled against phospho-ERK-1/2 (pT183/pY185) using immunofluorescence labeling with confocal microscopy imaging. (a, b) Immunofluorescence images are displayed together with (c) a surface rendering of the signal intensity within a suprabasal keratinocyte). The z-score normalized sampled signal intensity data (data points; plotted relative to the sample mean, \( \overline{\mu} \); and sample standard deviation, \( \overline{\sigma} \)) and loess smoothed signals (solid lines) associated with the (d) cytoplasm and (f) nuclei are displayed for Patient One (red), Two (green) and Three (blue), together with the 90 % confidence interval for positive and negative residuals (dashed lines), and the sampled data clouds for Patient One. Histograms of the discretized signal intensity within the (e) cytoplasm and (g) nuclei across individual spatial divisions are shown, together with the associated loess-curves (magenta line). The regions displayed in (b) and (c) are highlighted within (a) by the white and orange dashed lines, respectively. Scale bars represent 10 μm. Image data have undergone a non-linear transformation to improve visual appearance.

Additional file 4: Figure S5.

Human epidermis (Patient One) labeled against calmodulin using immunofluorescence labeling with confocal microscopy imaging. (a, b) Immunofluorescence images are displayed together with (c) a surface rendering of the signal intensity and (d) an isosurface volume rendering of a basal keratinocyte. The z-score normalized sampled signal intensity data (data points; plotted relative to the sample mean, \( \overline{\mu} \); and sample standard deviation, \( \overline{\sigma} \)) and loess smoothed signals (solid lines) associated with the (e) cytoplasm, (g) nuclei and (i) plasma-membrane are displayed for Patient One (red), Two (green) and Three (blue), together with the 90 % confidence interval for positive and negative residuals (dashed lines), and the sampled data clouds for Patient One. Histograms of the discretized signal intensity within the (f) cytoplasm, (h) nuclei and (j) plasma membrane across individual spatial divisions are shown, together with the associated loess-curves (magenta line). The regions displayed in (b) and (c, d) are highlighted within (a) by the white and orange dashed lines, respectively. Scale bars represent 10 μm. Image data have undergone a non-linear transformation to improve visual appearance.

Additional file 5: Figure S6.

Absolute thickness of epidermal tissue layers: (a) Basal Layer; (b) Spinous and Granular Layers; (c) Transitional Layer; (d) Total Epidermis. Some variation in the thickness of epidermal tissue layers was observed between patients. To examine this in a quantitative manner, a script was written in MATLAB to move along the lower boundary of each tissue layer and for every unique pixel to measure the minimum distance to the upper boundary (Fig. 2). These results were aggregated over all z-positions and target proteins/phospho-proteins, to estimate the distribution of epidermal tissue layer absolute thickness (μm) for Patient One (red), Two (green) and Three (blue). Note that the ‘peakiness’ of these distributions can be attributed to the aggregation of different z-positions and target proteins (with minor intra-patient variation for the epidermal thickness between tissue sections).

Additional file 6: Table S1.

Statistical associations between spatially-conditioned protein abundances used in this study. Statistical associations are italicized if they failed to exceed the data-derived significance thresholds of 0.469 for mutual information, and −0.284 and 0.286 for Pearson’s correlation. As shown in Fig. 3c & d, several relationships between the spatially-conditioned protein abundance data had a statistical association exceeded the data-derived threshold. The canonical interactions between (i) cytoplasmic phospho-Raf and phospho-MEK, and (ii) cytoplasmic phospho-MEK and phospho-ERK had relatively high mutual information and a positive Pearson’s correlation which was particularly strong for (ii). Relationships that reflect nucleocytoplasmic shuttling interactions were also relatively consistent with the known molecular translocation events with ERK-MAPK signal transduction; exceeding the data-derived thresholds for (iii) phospho-ERK-1/2 and (iv) phospho-MEK-1/2, but falling below these thresholds for (v) phospho-Raf-1.

Additional file 7: Figure S1.

A reaction kinetic scheme of the ERK-MAPK signaling cascade as modeled in this study. The spatial position within the epidermis modulates the relative abundance of Ca2+ and plasma-membrane calmodulin. Ca2+ activates Ras-GTP signaling, while CaM inhibits signal transduction from Ras-GTP to phosphorylated (pS338) Raf-1 (details given within the main text). The ERK-MAPK signaling cascade with Raf-1 phosphorylating and activating MEK-1/2, which phosphorylates ERK-1/2 is illustrated within the nucleus and cytoplasm; with nucleocytoplasmic shuttling reactions represented by corresponding reactions (dashed lines). Cytoplasmic phospho-ERK-1/2 phosphorylates Raf-1 to inhibit upstream signaling, while nuclear phospho-ERK-1/2 induces the expression of nuclear-localized DUSP4 which promotes ERK-1/2 dephosphorylation. Note that many of the species listed in this model were not measured, and due to the large number of unknown kinetic parameters associated with each reaction (including absolute phospho-protein concentrations), a simplified model was constructed using a normalized Hill differential equation approach (Fig. 1d).

Additional file 8: Figure S7.

Model simulation in the absence of calmodulin-mediated inhibition of Raf. To test the effects of removing calmodulin-mediated inhibition of phospho-Raf-1, the ‘Hill function weight parameter’ corresponding to this reaction was set to zero using SED-ML, and the model was evaluated at different spatial positions through the epidermis. Note that in comparison to Fig. 4, the relative abundance of phosphorylated ERK-MAPK is increased, and the profile of activation is much more linear, following the Ca2+ gradient (at top left). A SED-ML script to perform this simulation is included within the packaged SourceForge code.

Additional file 9: Figure S8.

Immunofluorescence data derived from whole-cell segmentation. The probability distribution function (p.d.f.) for phospho-MEK signal intensity within the (a) cytoplasm and (b) nucleus of segmented cells. Line color ranges from blue to red, corresponding to the increasing normalized distance value of nucleus centroids. Violin plots are also presented for the phospho-MEK signal intensity within the (c) cytoplasm and (d) nucleus of segmented cells, ordered along the x-axis by their normalized distance values (note the blue dashed vertical lines which demarcate the tissue layers, as labeled at bottom of (d), and green dashed lines which highlight ‘phospho-MEK bright basal cells’). Surface renderings of the p.d.f for phospho-MEK signal intensity within the cell (e) cytoplasm and (f) nucleus of segmented cells, plotted perpendicular to the normalized distance values of nuclei centroids (note the plots in surface renderings of the p.d.f are analogous to the ‘whole data clouds segmented by spatial position’, as shown in Additional file 2: Figure S3e and S3g in pMEK.png).

Additional file 10: Figure S10.

Gradient of keratinocyte cytoplasmic-to-nuclear ratio over the depth of the epidermis. A collection of interfollicular keratinocytes were segmented across multiple z-positions within several immunofluorescence data sets, as described in Materials and methods–Estimating cytoplasmic-to-nuclear volume ratio. A linear gradient of increasing cytoplasmic-to-nuclear volume ratio was applied across the depth of the epidermis, increasing from 2 to 5. Note the high degree of variance within the basal layer (Normalized Distance 0–1), which we believe is caused by proliferative cells growing at different stages of the cell cycle.

Additional file 11: Figure S9.

Human epidermis (Patient One) labeled against MEK-1/2 using immunofluorescence labeling with confocal microscopy imaging. (a, b) Immunofluorescence images are displayed together with (c) a surface rendering of the signal intensity within suprabasal keratinocytes. The z-score normalized sampled signal intensity data (data points; plotted relative to the sample mean, \( \overline{\mu} \); and sample standard deviation, \( \overline{\sigma} \)) and loess smoothed signals (solid lines) associated with the (d) cytoplasm and (f) nuclei are displayed for Patient One (red), Two (green) and Three (blue), together with the 90 % confidence interval for positive and negative residuals (dashed lines), and the sampled data clouds for Patient One. Histograms of the discretized signal intensity within the (e) cytoplasm and (g) nuclei across individual spatial divisions are shown, together with the associated loess-curves (magenta line). The regions displayed in (b) and (c) are highlighted within (a) by the white and orange dashed lines, respectively. Scale bars represent 10 μm. Image data have undergone a non-linear transformation to improve visual appearance. NB: the cytoplasmic fluorescence signal intensity (d) shows a strong decrease over the spinous and granular layers, particularly for Patient Two (green) and Three (blue). The inverse pattern for non-phosphorylated MEK-1/2 abundance relative to phospho-MEK-1/2 (Additional file 2: Figure S3d in pMEK.png) indicated that upon pS218/pS222 phosphorylation of MEK-1/2, there may have been masking of the epitope recognized by our total MEK-1/2 antibody. Testing this hypothesis is difficult and beyond the scope of our study, thus we excluded the non-phosphorylated/total ERK-MAPK components from the quantitative modeling.

Rights and permissions

Open Access  This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this licence, visit https://creativecommons.org/licenses/by/4.0/.

The Creative Commons Public Domain Dedication waiver (https://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Cursons, J., Gao, J., Hurley, D.G. et al. Regulation of ERK-MAPK signaling in human epidermis. BMC Syst Biol 9, 41 (2015). https://doi.org/10.1186/s12918-015-0187-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12918-015-0187-6

Keywords

BMC Systems Biology

ISSN: 1752-0509

Contact us