Employing a Mechanistic Model for the Mapk Pathway to Examine the Impact of Cellular all or None Behavior on Overall Tissue Response

Abstract

The mitogen activated protein kinase (MAPK) cascade is a three-tiered phosphorylation cascade that is ubiquitously expressed among eukaryotic cells. Its primary function is to propagate signals from cell surface receptors to various cytosolic and nuclear targets. Recent studies have demonstrated that the MAPK cascade exhibits an all-or-none response to graded stimuli. This study quantitatively investigates MAPK activation in Xenopus oocytes using both empirical and biologically-based mechanistic models. Empirical models can represent overall tissue MAPK activation in the oocytes. However, these models lack description of key biological processes and therefore give no insight into whether the cellular response occurs in a graded or all-or-none fashion. To examine the propagation of cellular MAPK all-or-none activation to overall tissue response, mechanistic models in conjunction with Monte Carlo simulations are employed. An adequate description of the dose response relationship of MAPK activation in Xenopus oocytes is achieved. Furthermore, application of these mechanistic models revealed that the initial receptor-ligand binding rate contributes to the cells' ability to exhibit an all-or-none MAPK activation response, while downstream activation parameters contribute more to the magnitude of activation. These mechanistic models enable us to identify key biological events which quantitatively impact the shape of the dose response curve, especially at low environmentally relevant doses.

1 INTRODUCTION

The mitogen-activated protein kinase (MAPK) cascades are an integral part of signal transduction in eukaryotic cells. These cascades are well-known, ubiquitous, and conserved in cells from yeast to mammals (Lewis et al., 1998; Widmann et al., 1999). In mammalian cells, there are three distinct, parallel cascades which have been well studied. These are the extracellular-signal-regulated-kinase (ERK) cascade, the c-Jun N-terminal kinase (JNK) cascade, and the p38 MAPK cascade. Each cascade is structured in a three-tier design, where the first level is a MAPK kinase kinase (MAPKKK), the second level is a MAPK kinase (MAPKK), and the final level is the MAPK (ERK, JNK or p38). A schematic of these cascades is depicted in Figure 1. Each kinase in the cascade phosphorylates and activates the immediate downstream kinase. The final component of the cascade, the MAPK, upon activation, phosphorylates various cytosolic targets and translocates to the nucleus to phosphorylate transcription factors affecting gene expression.

FIGURE 1:

Schematic of the three-tiered MAPK cascades (Adapted from Luster et al. (2004).

Studies have shown that the MAPK pathways can be activated by a wide variety of stimulants. Arsenic (Dong, 2002; Cooper et al., 2004), asbestos, cigarette smoke (Mossman et al., 2006), adenosine triphosphate (ATP) (Tai et al., 2001), copper (Mattie and Freedman, 2004), and numerous other compounds have been shown to be activators of the MAPK pathway. Similarly, activation of the MAPK pathway has been shown to lead to different end results. Transcriptional factors, such as AP-1, are activated (Shaulian and Karin, 2001; Zhong et al., 2005) and genes, such as c-fos and c-jun, are induced (Owuor and Kong, 2002). Additionally, depending on the amplitude and duration of the stimulus, MAPK activation can lead to cell proliferation or apoptosis (Ramos-Nino et al., 2002; Apati et al., 2003) in conjunction with other signaling cascades. Because the MAPK cascade can be activated by a variety of stimulants and plays such an important role in cell fate, it would be beneficial to mathematically model the signal transduction that occurs within the cascade.

Models of different aspects of the MAPK cascade have been presented throughout the literature. For example, El-Masri & Portier (1999) examined a mechanistic model of the MAPK pathway in relation to the cell cycle. It is well accepted that the MAPK pathway is governed by positive and negative feedback loops, and models incorporating the importance of this feedback have been described by Kholodenko (2007), Mutalik & Venkatesh (2006), and Salazar & Hofer (2006). These feedback loops have led to system models that exhibit bistability (Bhalla et al., 2002; Markevich et al., 2004; Ortega et al., 2006), ultrasensitivity (Kholodenko, 2000; Kolch et al., 2005), and oscillatory behavior (Chickarmane et al., 2007). A spatial model of the MAPK pathway was presented by Markevich et al. (2006), and Schwacke & Voit (2007) examined a model that accounted for interaction between parallel MAPK cascades. The aforementioned models are but a few examples; for a more detailed review of MAPK modeling efforts, see Orton et al. (2005) or Vayttanden et al. (2004).

Ferrell & Machleder (1998) and Mackeigan et al. (2005) postulate that MAPK phosphorylation occurs in an all-or-none, step-wise manner. From an initial glance at the dose response results exhibited by MAPK in the presence of an increasing stimulus, one might assume that the response is graded. That is, as the stimulus increases, the response increases uniformly. However, a graded tissue response could mean that each individual cell had a graded response, or it could be interpreted that each cell had a perfectly switch-like (all-or-none) response. Figure 2 displays a schematic for each of these possibilities. If the individual responses are graded, then each cell should have an increasing, intermediate level of MAPK phosphorylation as the stimulus increases. Under the switch like assumption, each cell can exhibit either a steady-state (none) level of MAPK phosphorylation, or a level which is elevated above the steady state level (all). Ferrell & Machleder (1998) present data that suggest that the latter, all-or-none response, is what occurs in the MAPK cascade. Studies suggest that the all-or-none behavior is also exhibited in gene induction and protein expression (Pirone and Elston, 2004; Zhang et al., 2006).

FIGURE 2:

Schematic of an all-or-none response (A) versus a graded response (B) in the presence of increasing stimulus.

Biologically based mathematical models are increasingly used to predict the shape of the dose response curve, specifically at low environmentally-relevant exposure levels. Hence, the ability to accurately model the all-or-none response becomes important in distinguishing between cell level and tissue level models. In this study, we investigate the all-or-none response using a mechanistic model that describes MAPK phosphorylation.

1.1 MATERIALS AND METHODS

All simulations were executed utilizing MATLAB® R2006b (The Mathworks, Inc., Natick, MA). The system of equations was integrated using MATLAB®'s ode15s command, a variable order method for solving a system of stiff differential equations. We examined MAPK phosphorylation mathematically using both empirical and mechanistic models. The empirical models utilized a Hill equation or a log normal cumulative distribution function. The mechanistic models were based on previously published models by El-Masri & Portier (1999) and Bhalla et al. (2002).

Each modeling effort was compared to MAPK phosphorylation data presented by Ferrell & Machleder (1998). The data were produced during an experiment in which Xenopus oocytes were exposed to varying amounts (ranging from 0 μM to 8 μM) of progesterone. The phosphorylation fold levels of MAPK were examined in 190 individual progesterone-treated oocytes, and 19 individual untreated oocytes. Each data point represents a sample of 11 to 39 oocytes. Data were digitized using the digitization software, TechDig v2.0 (Jones, R.B., Mudelein, Illinois, 1998).

All optimal parameters were identified using MATLAB®'s fminsearchbnd command (available at www.matlabcentral.com). This command is an implementation of the Nelder-Meade simplex algorithm which allows the user to specify lower and upper bounds for the possible parameter space. fminsearchbnd was used to find the parameters which minimize the following ordinary least squares (OLS) cost function (Sheiner and Beal, 1985),

J = \sum_{i = 1}^{N} | f (d_{i}; q) - z_{i} |^{2},

(1)

where f (d_i;q) is the model approximation at dose d_i, for the parameter set q, and zⁱ represents the data corresponding to dose d_i. Minimization of the cost function effectively minimizes the discrepancy between the model prediction and the available data.

2 EMPIRICAL MODELS FOR MAPK ACTIVATION

Ferrell & Machleder (1998) note that their data are approximated well by utilizing a Hill function with Hill coefficient nH = 1. The Hill function is given by the equation (Ferrell and Machleder, 1998)

C_{H i l l} (x) = \frac{x^{n H}}{k^{n H} + x^{a n H}}

(1)

where y represents the percentage of MAPK phosphorylation, x represents the amount of progesterone, nH is the Hill coefficient, and k represents the concentration of x at which the oocytes respond half-maximally. This Hill equation fits the data with a least squares cost function of 1519, and its approximation to the data is depicted in Figure 3. This Hill equation approximates MAPK activation at a tissue level, and thus, displays a graded response.

FIGURE 3:

Comparison of a hill equation and a lognormal cumulative distribution function approximation of MAPK phosphorylation.

Under the all-or-none MAPK activation assumption, the level of stimulus required to fully phosphorylate each cell varies (i.e. one cell may exhibit full MAPK phosphorylation at a different dose than another). We further investigated the impact of the distribution of activation doses on the dose response curve. Studies have shown that the uptake of radioactivity is log-normally distributed among a group of cells (Neti and Howell, 2006). If we consider a similar assumption for MAPK activation, that is, that the amount of stimulus required to cause full phosphorylation is distributed log-normally among the cells, then the percentage of cells which have become activated for any given dose can be represented by the cumulative distribution function for a lognormal distribution. This cumulative distribution function takes the form

C_{L o n g - n o r m} (x) = \frac{1}{2} + \frac{1}{2} e r f [\frac{ln (x) - μ}{σ ✓ 2}],

(2)

where x is the level of stimulus, μ is the logarithmic mean of the activation dose, σ is the logarithmic standard deviation, and erf is the error function, given by

\begin{matrix} e r f (x) = & \frac{2}{✓ π} & \int_{0}^{x} e^{- t^{2}} d t \end{matrix}

(3)

Using the fminsearch optimization routine, we identify optimal values for the parameters μ and σ. The best fit to the data was produced when μ = −3.43 (μM) and σ = 2.78 (μM). Note that since we are looking at the log normal mean and standard deviation, the optimized value of μ does not reflect a negative value, but rather a value on the order of 0.001. This approximation results in a least squares cost function of 821, and is shown in comparison to the Hill approximation and the data in Figure 3. Similar to the Hill equation approximation, this lognormal cumulative distribution function represents MAPK activation at the tissue level and exhibits a graded response.

While both the Hill equation and the lognormal cumulative distribution function give reasonable fits to the Ferrell & Machleder (1998) data, they both are representative of a tissue level graded response. However, experimental evidence has shown that on the cellular level an all-or-none response is exhibited (Ferrell and Machleder, 1998; Bagowski and Ferrell, 2001). To adequately investigate the role of MAPK in Xenopus oocytes, it would be beneficial to have the ability to focus on both tissue level and cellular level responses. The empirical models presented here lack description of key biological processes at the cellular level and therefore give no insight into whether the cellular response occurs in a graded or all-or-none fashion. In this paper, we employ mechanistically-based models of MAPK phosphorylation to investigate the propagation from cellular to tissue response.

3 EMPLOYING A MECHANISTIC MAPK MODEL

3.1 Model Structure

We consider a model that accounts for the biological structure of the MAPK signaling cascade. The mechanistic model which we initially consider is based upon a model of the MAPK cascade developed by El-Masri & Portier (1999). This model, whose end result predicts fold levels of AP-1 production, consists of three sections. The first section, dealing with the activation of protein kinase C (PKC), focuses on the initial interpretation of the external signal. The second section, which involves the activation of MAPK kinases, models the translation of the growth signal from the cytoplasm to DNA. The final section deals with the nuclear events leading to the expression of AP-1. The primary objective of the study by El-Masri and Portier (1999) was to model the effect of MAPK on cellular replication, hence the inclusion of the transcripton factor AP-1. Because the main focus of our study deals with MAPK, only the first two sections of the model will be utilized. In the following, we give a brief detailed description of the activities occurring during each section of the model. For a more in-depth description, we refer readers to El-Masri & Portier (1999). Model equations are presented in Appendix A, and parameter values are reported in Tables I and II.

The first section of the model has been adapted from a previously published calcium oscillator model (Cuthbertson and Chay, 1991). The signaling cascade is activated when growth factors bind to cellular receptors. This binding initiates a build up of activated guanosine triphosphate (GTP) binding proteins, which, in turn, leads to the activation of phospholipid lipase C (PLC). This causes the hydrolysis of phosphatidylinositol (4,5)-bisphosphate (PIP₂), producing inositol (1,4,5) triphosphate (InsP₃) and diacylglycerol (DAG). InsP₃ facilitates the release of Ca₂+ from the endoplasmic reticulum (ER). This free Ca₂+, along with DAG, activates PKC. Activated PKC is then responsible for degrading active GTP (a negative feedback process), as well as initiating the next portion of the model.

In the second section of the model, activation of PKC is transferred into the nucleus by multiple steps involving phosphorylation and dephosphorylation of MAPK (Yamaguchi et al., 1995). Initially, PKC activates Raf, which then activates MAP kinase kinase (MAPKK). Activation of MAPKK ultimately increases nuclear MAPK activity. A negative feedback process is present, as MAPK deactivates Raf, essentially regulating the activation of the cascade.

The model is used to generate time course data, as well as a dose response curve. Time course and dose response results generated using the parameter values listed in Tables I and II are displayed in Figures 4 and 5, respectively. The dose response results are illustrated by the solid blue line in Figure 5. The system was studied under different stimulus levels to examine the presence of an all-or-none activity. Figure 4 was produced by running the model with a constant stimulus (of 50 nM or 800 nM) present from time t = 1000 min to t = 1200 min. The first 1000 simulated minutes were run stimulus-free in order to attain steady state response folds of unity prior to the introduction of the stimulus. Note that a dose of 50 nM leads to negligible activation above background (less than 1.01 fold induction), while the fold level of phosphorylated MAPK elevates in the presence of the 800 nM stimulus. We believe that the 50 nM level does not elicit an elevation of phosphorylated MAPK because there is not enough stimulus present to overcome the negative feedback built into the model. This negative feedback, inherent to the system, keeps the low dose exposure response close to the steady state levels.

FIGURE 4:

Time course results for MAPK fold induction, with stimulus present from t=1000 to t=1200 minutes. Note that the model was run for 1000 simulated minutes (stimulus free) to establish a steady state.

FIGURE 5:

Sensitivity comparison for the parameter V_mapk, displayed as MAPK fold induction at t=1100 min. The solid blue line represents model results using the original parameter value.

TABLE I:

Parameter values (Cuthbertson and Chay, 1991) for the PKC activation portion of the mechanistic model

Parameter	Function	Value	Units
K_GF	1^st order constant of stimulus binding to receptor	1	min⁻¹
hg	rate constant of GTP inactivation	300	min⁻¹
K_d	Hill binding coefficient of DAG to PLC	5	nM
K_g	Hill binding coefficient of GTP to PLC	50	nM
n	Hill exponent – DAG binding to PLC	2	unitless
m	Hill exponent – GTP binding to PLC	4	unitless
K_DAG	rate of DAG formation by PLC	30000	nM/min
hd	rate constant for breakdown of DAG	240	min⁻¹
ld	rate of background release of DAG	18	nM/min
K_s	Hill binding coefficient of InsP₃ to calcium membrane	5	nM
hc	rate constant of calcium uptake	24	min⁻¹
lc	background production rate of calcium	4800	nM/min
K_P	Binding coefficient for DAG to PKC	40	nM
K_C	Binding coefficient for calcium to PKC	100	nM

TABLE II:

Parameter values (El-Masri and Portier, 1999) for the MAPK activation section of the mechanistic model. Values denoted with a (*) have been modified from the original sources.

Parameter	Function	Value	Units
RPKCⁱ	background activation of PKC	0.006	unitless
V_raf	Phosphorylation rate constant of Raf	1.6	min⁻¹
K_dc	rate constant for deactivation of Raf via RMAPK	0.04	min⁻¹
V_mapkk	rate constant for phosphorylation of MAPKK	0.02	min⁻¹
V_dmapkk	rate constant for MAPKK* deactivation	2	min⁻¹
V_mapk	rate constant for phosphorylation of MAPK	0.01	min⁻¹
V_dmapk	rate constant for MAPK* deactivation	1.1	min⁻¹
α₁	ratio used to determine Raf steady state	0.2247*	unitless
α₂	ratio used to determine MAPKK steady state	0.0099*	unitless
α₃	ratio used to determine MAPK steady state	0.009*	unitless

The dose response data given in Figure 5 was generated by running the model with varying stimulus levels, and plotting the fold induction level for a given time (t = 1100) at each level. As the stimulus level increases, the fold induction level reaches a higher steady state. Notice that fold levels are close to the unstimulated steady state at doses less than 100 nM. Again, this may be attributed to the negative feedback in the system.

In order to better understand the role that the parameters play in the model, a visual sensitivity analysis was conducted. Parameter values were varied by a factor of 2, and the ensuing changes in the dose response curve of phosphorylated MAPK was examined. Figures 5 and 6 display the sensitivity analysis results for a few selected parameters. V_mapk (in Figure 5) represents the rate constant for MAPK activation. Examination of this figure shows that the parameter (analogous to other activation parameters) has a marked effect on the maximum fold induction level of phosphorylated MAPK reached. The results shown in Figure 5 are representative of the results for the other activation parameters (V_mapkk, V_raf, etc.). The parameter whose results are shown in Figure 6, K_GF, represents the first order rate constant which governs the rate at which the stimulus binds to the receptor, thus activating the cascade. We notice that changes in the value of K_GF shifts the dose response curve horizontally. These horizontal shifts impact the lowest levels of stimulus at which MAPK activation increases over background. It should be noted that this is the only parameter in the system which exhibits such an effect in the sensitivity analyses. Thus, variations of the value of K_GF reflect an increased or decreased sensitivity of the cell to the stimulus.

FIGURE 6:

Sensitivity analysis for K_GF, displayed as MAPK fold induction at t=1100 min. K_GF = 1 represents the model results using the original parameter value.

In this study, we consider the percent activation of MAPK in tissue. Our mechanistic model predicts fold induction of activated MAPK in cells, which is converted to percent activation. Then, the steady state level of unity will correspond to zero percent activation, while the maximum fold induction corresponds to 100% activation. This will normalize the fold induction results. Thus, the variation in the maximum amplitude attained in Figure 5 would be rendered insignificant, as each maximum would correspond to 100% activation. Therefore, the various activation parameters are of little concern to us, allowing us to concentrate our efforts on the K_GF parameter.

3.2 Monte Carlo Simulation

The results produced by the mechanistic model can be considered the MAPK fold level for a single cell, or it can represent the tissue-level response exhibited under the assumption that each cell demonstrates a graded response. In order to implement this model with the all-or-none assumption (with single cells either responding or not responding), we assume that the model results are representative of a single cell, and we employ a Monte Carlo simulation, as follows.

We first assume that the K_GF parameter varies from cell-to-cell and is log-normally distributed (with mean μ _K _GF and standard deviation σ _K _GF ). Assigning variability to the parameter K_GF can be supported by several mechanisms. First, it is possible that the number of receptors varies from cell to cell, causing some cells to be more sensitive than others. Perhaps there is a structural difference in the receptors, allowing some receptors to have a higher affinity to the stimulus. Finally, there could be a difference in the concentration of stimulus reaching the cells. While these factors are speculative, they do provide motivation for varying K_GF in the model.

The model is evaluated for a sampling of K_GF in the distribution. If the fold level of phosphorylated MAPK is elevated (as illustrated by the example of the 800 nM results in Figure 4), then the cell is considered to be on. If the fold level remains equal to 1, then the cell is off. The percent of activated MAPK is defined to be the percentage of activated cells (number of active cells ÷ total number of cells × 100).

The Monte Carlo simulation of the MAPK model was run with an arbitrary number of n=100 and optimized to find the logarithmic mean and standard deviation for K_GF, using the data from Ferrell & Machleder (1998). We found optimal parameters μ _K _GF = 0.77 (1/min) and σ _K _GF = 2.35 (1/min), resulting in a least squares cost function value of 832. The resulting MAPK activation is illustrated by the solid blue line in Figure 7 in comparison to the data.

FIGURE 7:

Approximation of MAPK activation computed using the Monte Carlo simulation with unbounded (solid blue line) and bounded (dashed green line), log normally distributed parameter KGF compared to data reported by Ferrell & Machleder (1998).

Upon examination of these results in Figure 7, we notice that the low dose (< 1 nM) results are not in accord with our expectations. As evidenced by the data, we would expect that low doses would have little effect on MAPK activation. These initial results, however, suggest that at very low doses, 8% of the cells exhibit MAPK activation. By assuming that the parameter K_GF has a log normal distribution, we subject our model to outliers which have the potential to generate a few values that are much larger or much smaller than the majority of the others. We believe that this is the cause for the discrepancy between simulated results and data at low doses.

If we instead consider a log normal distribution (using the optimal parameters μ _K _GF = 0.77 (1/min) and σ _K _GF = 2.35 (1/min)) that is bounded above and below (i.e. K_GF ∊ [0.05, 40]), then our model produces the results depicted by the dashed green line in Figure 7. The corresponding cost function value is now J = 701, implying that the fit is improved, especially when considering the low dose data. The need to limit the possible values of K_GF is not surprising, as physiological parameters need lie in a range that is biologically feasible.

3.3 Further Investigation

Bhalla et al. (2002) present a more complex model of the MAPK cascade. Their model is comprised of 105 binding and enzyme reactions, modeled by a system of 90 differential equations. The model was developed to examine the role of MAP kinase phosphatase on the MAPK cascade. Phosphorylation of the MAPKs occurs as a dual step process, and forward and backward feedback loops are inherent in the model.

While there are many differences between the Bhalla et al. (2002) model and our MAPK model, there are also several similarities. An end result of both models is the prediction of phosphorylated MAPK. Also, activation of the cascade is initiated by a binding process in each model. That is, there is a binding coefficient in the Bhalla et al. (2002) model that is analogous to the K_GF parameter of our model. A brief sensitivity analysis was conducted, varying the value of the initial binding parameter, and the results are presented in Figure 8. These results appear to be similar to those presented in Figure 6 for the El-Masri and Portier (1999) Model.

FIGURE 8:

Sensitivity analysis for binding parameter in the Bhalla et al. (2002) model.

To further examine our method for incorporating the all-or-none response into a mechanistic model, we employ the previously described Monte Carlo method in conjunction with the Bhalla et al. (2002) model. We assume that the binding coefficient from the Bhalla et al. (2002) model has the same distribution that we found to be optimal for our model. Thus, the binding coefficient is log-normally distributed with a mean of 0.77 (1/min) and a standard deviation of 2.35 (1/min), and is contained in the interval [0.05, 40]. The Bhalla et al. (2002) Model is evaluated with each value in a sampling (n = 100) of the distribution. If the level of phosphorylated MAPK increases above background, then the cell is considered to be “on”. The percent of activated MAPK is again defined to be the percentage of activated cells. Results generated by the Bhalla et al. (2002) model are shown in comparison to the results generated by our model, as well as the data from Ferrell & Machleder (1998), in Figure 9. The least squares cost function value associated with the Bhalla et al. model is 869. Examining Figure 9, we see that the use of our optimal distribution parameters with the Bhalla et al. (2002) model produces similar results (which could be improved by finding optimal parameters for it, rather than using parameters optimized for the El-Masri & Portier (1999) model). This suggests that our Monte Carlo method for incorporating the all-or-none response can be employed in conjunction with different mechanistic models.

FIGURE 9:

Approximation of MAPK activation under the assumption of a cellular all or none response, computed using the Monte Carlo simulation in conjunction with the Bhalla et al. (2002). Model compared to results from the El-Masri model (El-Masri and Portier, 1999) and data reported by Ferrell & Machleder (1998).

A comparison of the least squares cost function values, presented in Table III, shows that the mechanistic models, modified to incorporate the all-or-none MAPK response, provide a better fit to the data than the empirical models do. The benefit of using these mechanistic models is that they can incorporate more information about the biological processes. They also help to provide an understanding of key biological events contributing to the shape of the dose response curve. In this example, we found that the initial receptor-ligand binding rate contributes to the cell's ability to exhibit an all-or-none MAPK activation response, while downstream activation parameters contribute more to the magnitude of activation (see Figure 5).

TABLE III:

Optimal parameters and corresponding cost values

Model Description	Optimal Parameters	OLS Cost Value
Hill Model^a	K = 31.25	1519
Lognormal distribution	μ = −3.43, σ = 2.78	821
El-Masri, All-or-None, unbounded^a	μ _KGF = 0.77, σ _KGF = 2.35	832
El-Masri, All-or-None, bounded^b	μ _KGF = 0.77, σ _KGF = 2.35, K_GF ∊ [0.05 40]	701
Bhalla, All-or-None, bounded^{c *}	μ _K = 0.77, σ _K = 2.35, K ∊ [0.05 40]	869
El-Masri, Graded^b	K_GF = 23.98	2580
Bhalla, Graded^c	K = 5.66	3127

model from Ferrell and Machleder (1999);

model from El-Masri and Portier (1999);

model from Bhalla et al. (2002);

parameters for the Bhalla, All-or-None model were not optimized

3.4 All-or-None Cellular Response vs. Graded Cellular Response

When coupled with the Monte Carlo simulation, the El-Masri & Portier (1999) and Bhalla et al. (2002) MAPK models can predict MAPK activation, assuming an all-or-none cellular response. These models, without the Monte Carlo simulation, can also be used to predict overall MAPK activation if we assume that the cellular response is graded. Under the assumptions that each cell exhibits an intermediate graded response to increasing stimulus and that the tissue level response is an average of the cellular responses, then the results generated by these models (without the Monte Carlo simulation) would represent tissue level MAPK activation with graded cellular responses.

Again using the data reported in Ferrell & Machleder (1998), we wish to compare the accuracy of the graded response model with our implementation of the all-or-none model. Since the El-Masri & Portier (1999) and Bhalla et al. (2002) models predict fold induction and actual units of MAPK activation, respectively, the results of both models are linearly scaled so that the maximum induction relates to a 100% activation level and basal levels become 0% activation. Results generated using the models with a graded cellular response are shown in Figure 10. Optimization of the graded models to the Ferrell & Machleder (1998) data yields a binding parameter of 5.66 sec⁻¹ for the Bhalla et al. (2002) model and a K_GF value of 23.98 min⁻¹ for the El-Masri & Portier (1999) model. The optimal parameter values produce cost function values of 3127 and 2580, respectively. These cost function values are considerably higher than those computed using the models with an all-or-none approach, suggesting that the ability of the models to predict MAPK activation is enhanced by assuming that cellular MAPK responses occur in an all-or-none fashion.

FIGURE 10:

MAPK activation results predicted by the Bhalla et al. (2002) Model and the El-Masri & Portier (1999) model for a graded cellular response compared to data reported by Ferrell & Machleder (1998).

4. DISCUSSION

Many signaling biological processes (nuclear transcriptional receptors, kinase or phosphatase cascades, G-coupled protein receptors, etc.) are mediated through cascading processes leading to cellular events (e.g., commitment to cell division, cell differentiation, and phenotypic alterations) that eventually translate to tissue and/or systemic responses. These processes have been typically examined and modeled at a tissue level, but technology has facilitated the examination of these events at the cellular level.

Experimental evidence of either a tissue-based graded or cellular based all-or-none response was investigated using Xenopus oocyte maturation in the presence of progesterone (Ferrell and Machleder, 1998). For populations of Xenopus oocytes, the progesterone concentrations leading to maturation were distributed fairly broadly. Data for the tissue response were modeled empirically, and it was found that the data were fit well by a Hill equation with Hill coefficient of approximately 1, implying a graded response. However, Hill coefficients for individual oocytes were found to be between 20 and 30 (Ferrell and Machleder, 1998). These high values indicate an all-or-none switch for maturation in any individual oocyte.

While several models of varying levels of complexity have been developed for the MAPK signaling cascade (see Orton et al. (2005) for a detailed review), few, if any, of the published models account for the all-or-none response. In this paper, we utilized biological knowledge of underlying mechanisms for individual oocyte maturation using a detailed description of the MAPK pathway to investigate the impact of all-or-none behavior of individual cells on the overall oocyte population dose-response relationship. We have presented a Monte Carlo simulation approach which can be used to incorporate the all-or-none behavior into MAPK modeling efforts. Using this approach with our MAPK model (El-Masri and Portier, 1999), we recapitulated experimental data presented by Ferrell & Machleder (1998). We further investigated our method with a more complicated MAPK model developed by Bhalla et al. (2002). Both models, in conjunction with the Monte Carlo simulation, generated results which are consistent with the Ferrell & Machleder (1998) data. Furthermore, in both models, it was determined that the critical parameter for incorporating the all-or-none response was the initial binding rate parameter. Additionally, we compared model results generated with the assumption that cellular responses occur in a graded manner with results compiled under the assumption that an all-or-none cellular response is present. We found that the model incorporating the all-or-none response provides a better fit to MAPK activation data, then our implementation of the graded response model.

The examined importance of the initial binding rate parameter suggested that receptor binding was the driving mechanism for the presence of a threshold-like level of the ligand to increase the activity of MAPK leading to oocyte maturation. As noted earlier, allowing the binding rate to vary across cells can potentially be explained by several factors including availability of the ligand, affinity of the ligand to the receptors, structure of the receptor and quantity of available receptor sites. Our model-driven investigation provides an example where mechanistic information relating to the MAPK pathway's role in cell maturation was organized based on biological representation of important processes to investigate dose-response at the cellular level. This mechanistic description of the cellular-response can be introduced into a larger computational model where receptor binding via the MAPK pathway is part of a general mode of action for toxicological effects (such as cellular proliferation) to describe an overall biologically based dose-response (BBDR) model.

Examples of signal-mediated cascading biological motifs are plentiful. These signaling pathways have been described experimentally and computationally either as “all-or-none” phenomena or as graded responses. For example, neuronal cell firing in response to neurotransmitters is an all-or-none phenomenon, where as biochemical changes in response to activation of the steroid receptor superfamily, such as the progesterone receptor, have been described as graded responses. As technology has advanced, many of the observed graded responses appear due to an underlying all or none response at the cellular level. For example, the transcriptional changes due to activation of hepatic nuclear receptors, such as the Ah receptor or the Constitutive Androstane Receptor, once thought to yield graded responses at the tissue level, are due to dose dependent changes in the number of fully activated cells (Buhler et al., 1992; Tritscher et al., 1992; Fox et al., 1993). In many cases, these cascading biological motifs have nonlinear dose response behaviors for endogenous ligands and for exogenous toxicants, acting as switches with “all-or-none” responses over a narrow range of concentrations (Andersen et al., 2002). In addition, these induction responses appear to occur at the level of individual cells, thus a 50% response of the liver for induction does not represent 50% induction in each individual cell. Instead, half of the cells are fully induced and half are unaffected, while cells “switch” from one phenotypic state to another (Andersen, 2002). The difference between a graded response and a switch-like behavior may have significant implications for understanding dose response relationships for endogenous hormones or exogenous toxicants.

Biologically based dose response models are used to provide simulations that link molecular events to tissue and organ level responses such as morphogenesis and carcinogenesis. These models, such as the one we used for MAPK activation, accounting for the biological switches, would improve our understanding of these important biological processes. Similarly, BBDR models may prove useful in linking molecular, cellular, tissue and organ responses which will predict adverse health effects. These BBDR models must account for normal control of the signaling motifs and for the perturbations by toxic compounds (Hoel et al., 1983; Andersen, 2002). As the biological basis of these switches becomes unraveled and incorporated into computational models, these BBDR models should eventually serve to improve risk assessment with a variety of toxicants exhibiting receptor-based modes of action.

Footnotes

ACKNOWLEDGMENTS

This manuscript has been reviewed in accordance with the policy of the National Health and Environmental Effects Research Laboratory, U.S. Environmental Protection Agency, and approved for publication. Approval does not signify that the contents necessarily reflect the views and policies of the Agency, nor does mention of trade names or commercial products constitute endorsement or recommendation for use.

The authors would like to thank Dr. Karen Yokley for valuable insight provided during the completion of this project. Additionally, the authors express appreciation towards Dr. Marina Evans, Dr. Rory Conolly, and Dr. Linda Birnbaum for helpful comments during the preparation of this manuscript.

APPENDIX A

References

Andersen

M. E.

. (2002). The use of quantitative histological and molecular data for risk assessment and biologically based model development. Toxicol Pathol 30, 106–111.

Andersen

M. E.

Yang

R. S.

French

C. T.

Chubb

L. S.

, and Dennison

J. E

. (2002). Molecular circuits, biological switches, and nonlinear dose-response relationships. Environ Health Perspect 110 Suppl 6, 971–978.

Apati

Janossy

Brozik

Bauer

P. I.

, and Magocsi

. (2003). Calcium induces cell survival and proliferation through the activation of the MAPK pathway in a human hormone-dependent leukemia cell line, TF-1. J Biol Chem 278, 9235–9243.

Bagowski

C. P.

, and Ferrell

J. E.

Jr. (2001). Bistability in the JNK cascade. Curr Biol 11, 1176–1182.

Bhalla

U. S.

Ram

P. T.

, and Iyengar

. (2002). MAP kinase phosphatase as a locus of flexibility in a mitogen-activated protein kinase signaling network. Science 297, 1018–1023.

Buhler

Lindros

K. O.

Nordling

Johansson

, and Ingelman-Sundberg

. (1992). Zonation of cytochrome P450 isozyme expression and induction in rat liver. Eur J Biochem 204, 407–412.

Chickarmane

Kholodenko

B. N.

, and Sauro

H. M

. (2007). Oscillatory dynamics arising from competitive inhibition and multisite phosphorylation. J Theor Biol 244, 68–76.

Cooper

K. L.

Myers

T. A.

Rosenberg

Chavez

, and Hudson

L. G

. (2004). Roles of mitogen activated protein kinases and EGF receptor in arsenite-stimulated matrix metalloproteinase-9 production. Toxicol Appl Pharmacol 200, 177–185.

Cuthbertson

K. S.

, and Chay

T. R

. (1991). Modelling receptor-controlled intracellular calcium oscillators. Cell Calcium 12, 97–109.

10.

Dong

. (2002). The molecular mechanisms of arsenic-induced cell transformation and apoptosis. Environ Health Perspect 110 Suppl 5, 757–759.

11.

El-Masri

H. A.

, and Portier

C. J

. (1999). Replication potential of cells via the protein kinase C-MAPK pathway: application of a mathematical model. Bull Math Biol 61, 379–398.

12.

Ferrell

J. E.

Jr. , and Machleder

E. M

. (1998). The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes. Science 280, 895–898.

13.

Fox

T. R.

Best

L. L.

Goldsworthy

S. M.

Mills

J. J.

, and Goldsworthy

T. L

. (1993). Gene expression and cell proliferation in rat liver after 2,3,7,8-tetrachlorodibenzo-p-dioxin exposure. Cancer Res 53, 2265–2271.

14.

Hoel

D. G.

Kaplan

N. L.

, and Anderson

M. W

. (1983). Implication of nonlinear kinetics on risk estimation in carcinogenesis. Science 219, 1032–1037.

15.

Kholodenko

B. N

. (2000). Negative feedback and ultrasensitivity can bring about oscillations in the mitogen-activated protein kinase cascades. Eur J Biochem 267, 1583–1588.

16.

Kholodenko

B. N

. (2007). Untangling the signalling wires. Nat Cell Biol 9, 247–249.

17.

Kolch

Calder

, and Gilbert

. (2005). When kinases meet mathematics: the systems biology of MAPK signalling. FEBS Lett 579, 1891–1895.

18.

Lewis

T. S.

Shapiro

P. S.

, and Ahn

N. G

. (1998). Signal transduction through MAP kinase cascades. Adv Cancer Res 74, 49–139.

19.

Mackeigan

J. P.

Murphy

L. O.

Dimitri

C. A.

, and Blenis

. (2005). Graded mitogen-activated protein kinase activity precedes switch-like c-Fos induction in mammalian cells. Mol Cell Biol 25, 4676–4682.

20.

Markevich

N. I.

Hoek

J. B.

, and Kholodenko

B. N

. (2004). Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades. J Cell Biol 164, 353–359.

21.

Markevich

N. I.

Tsyganov

M. A.

Hoek

J. B.

, and Kholodenko

B. N

. (2006). Long-range signaling by phosphoprotein waves arising from bistability in protein kinase cascades. Mol Syst Biol 2, 61.

22.

Mattie

M. D.

, and Freedman

J. H

. (2004). Copper-inducible transcription: regulation by metal- and oxidative stress-responsive pathways. Am J Physiol Cell Physiol 286, C293–301.

23.

Mossman

B. T.

Lounsbury

K. M.

, and Reddy

S. P

. (2006). Oxidants and signaling by mitogen-activated protein kinases in lung epithelium. Am J Respir Cell Mol Biol 34, 666–669.

24.

Mutalik

V. K.

, and Venkatesh

K. V

. (2006). Effect of the MAPK cascade structure, nuclear translocation and regulation of transcription factors on gene expression. Biosystems 85, 144–157.

25.

Neti

P. V.

, and Howell

R. W

. (2006). Log normal distribution of cellular uptake of radioactivity: implications for biologic responses to radiopharmaceuticals. J Nucl Med 47, 1049–1058.

26.

Ortega

Garces

J. L.

Mas

Kholodenko

B. N.

, and Cascante

. (2006). Bistability from double phosphorylation in signal transduction. Kinetic and structural requirements. FEBS J 273, 3915–3926.

27.

Orton

R. J.

Sturm

O. E.

Vyshemirsky

Calder

Gilbert

D. R.

, and Kolch

. (2005). Computational modelling of the receptor-tyrosine-kinase-activated MAPK pathway. Biochem J 392, 249–261.

28.

Owuor

E. D.

, and Kong

A. N

. (2002). Antioxidants and oxidants regulated signal transduction pathways. Biochem Pharmacol 64, 765–770.

29.

Pirone

J. R.

, and Elston

T. C.

. (2004). Fluctuations in transcription factor binding can explain the graded and binary responses observed in inducible gene expression. J Theor Biol 226, 111–121.

30.

Ramos-Nino

M. E.

Haegens

Shukla

, and Mossman

B. T

. (2002). Role of mitogen-activated protein kinases (MAPK) in cell injury and proliferation by environmental particulates. Mol Cell Biochem 234–235, 111–118.

31.

Salazar

, and Hofer

. (2006). Kinetic models of phosphorylation cycles: a systematic approach using the rapid-equilibrium approximation for protein-protein interactions. Biosystems 83, 195–206.

32.

Schwacke

J. H.

, and Voit

E. O

. (2007). The potential for signal integration and processing in interacting MAP kinase cascades. J Theor Biol 246, 604–620.

33.

Shaulian

, and Karin

. (2001). AP-1 in cell proliferation and survival. Oncogene 20, 2390–2400.

34.

Sheiner

, and Beal

. 1985. Pharmacokinetic parameter estimates from several least squares procedures: Superiority of extended least squares. J Pharmacokinet Pharmacodyn 13(2), 185–201.

35.

Tai

C. J.

Kang

S. K.

Tzeng

C. R.

, and Leung

P. C

. (2001). Adenosine triphosphate activates mitogen-activated protein kinase in human granulosa-luteal cells. Endocrinology 142, 1554–1560.

36.

Tritscher

A. M.

Goldstein

J. A.

Portier

C. J.

McCoy

Clark

G. C.

, and Lucier

G. W

. (1992). Dose-response relationships for chronic exposure to 2,3,7,8-tetrachlorodibenzo-p-dioxin in a rat tumor promotion model: quantification and immunolocalization of CYP1A1 and CYP1A2 in the liver. Cancer Res 52, 3436–3442.

37.

Vayttaden

S. J.

Ajay

S. M.

, and Bhalla

U. S

. (2004). A spectrum of models of signaling pathways. Chembiochem 5, 1365–1374.

38.

Widmann

Gibson

Jarpe

M. B.

, and Johnson

G. L

. (1999). Mitogen-activated protein kinase: conservation of a three-kinase module from yeast to human. Physiol Rev 79, 143–180.

39.

Yamaguchi

Ogita

Nakamura

, and Nishizuka

. (1995). The protein kinase C isoforms leading to MAP-kinase activation in CHO cells. Biochem Biophys Res Commun 210, 639–647.

40.

Zhang

Andersen

M. E.

, and Conolly

R. B.

(2006). Binary gene induction and protein expression in individual cells. Theor Biol Med Model 3, 18.

41.

Zhong

C. Y.

Zhou

Y. M.

Douglas

G. C.

Witschi

, and Pinkerton

K. E

. (2005). MAPK/AP-1 signal pathway in tobacco smoke-induced cell proliferation and squamous metaplasia in the lungs of rats. Carcinogenesis 26, 2187–2195.