# Estimating Network Kinetics of the MAPK/ERK Pathway Using Biochemical Data

Authors: Vilda Purutçuoğlu; Ernst Wit
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

The MAPK/ERK pathway is a major signal transduction system
which regulates many fundamental cellular processes including the
growth control and the cell death. As a result of these roles, it has
a crucial importance in cancer as well as normal developmental processes.
Therefore, it has been intensively studied resulting in a wealth
of knowledge about its activation. It is also well documented that
the activation kinetics of the pathway is crucial to determine the nature
of the biological response. However, while individual biochemical
steps are well characterized, it is still difficult to predict or even understand
how the activation kinetics works. The aim of this paper is to
estimate the stochastic rate constants of the MAPK/ERK network dynamics.
Accordingly, taking a Bayesian approach, we combined underlying
qualitative biological knowledge in several competing dynamic
models via sets of quasireactions and estimated the stochastic rate constants
of these reactions. Comparing the resulting estimates via the
BIC and DIC criteria, we chose a biological model which includes
EGFR degradation—Raf-MEK-ERK cascade without the involvement
of RKIPs.


---

## Body

## 1. Introduction

All cellular responses are regulated by various signal transduction pathways. The signal transduction starts by an external stimulus, usually a ligand binding to a receptor at the cell surface. They generate intracellular signals that are transmitted and integrated through biochemical reactions. These biochemical reactions often include changes in gene expressions and particular biological responses, which can include the cell reproduction, the motility, and others. On the other hand, any malfunction in these mechanisms has a direct influence on the expression or on the function of gene products which are components of these regulatory mechanisms, hereby resulting in alterations of the biological responses and many illnesses such as heart disease, diabetes, and cancer [1]. Hence, knowledge about these pathways is very helpful for understanding the behaviour of biological activities and identifying drug targets which can be seen as the major aim of most of the biochemical and bioengineering studies.The MAPK (mitogen-activated protein kinase) also known as ERK (extracellular signal-regulated kinase) pathway is one of the major signal transduction systems which are involved in the cellular growth control of all eukaryotes from the cell reproduction to death. The basic structure of the MAPK/ERK pathway includes a number of phosphorylations on the protein level. These phosphorylations are directed by core components which are Ras, Raf, MEK, and ERK proteins and regulatory components, such as ERK, RKIP, and RSK proteins. The phosphorylations by core components transduce the signal from the cell membrane to nucleus and the phosphorylations via regulatory components modulate the efficiency and the duration of the signal transduction through the pathway. The effects of these components are conducted by positive and negative feedback loops. The functionality of the MAPK pathway gives rise to nonlinear behaviour, such as ultrasensitivities, bistabilities, and periodic behaviour. Moreover, there is an influence of the effect of scaffolding proteins and subcellular compartmentalisation which makes the functionality of the system even more complex [2]. Because of these underlying features, besides the execution of the large number of interactions at the protein level, the outcomes of the signal transfer are stochastic in nature. Thus, the interaction maps of the proteins or the simple representation of the system like given in Figure 4 at Appendix are not enough to make predictive statements about this network structure, although they do provide us with a starting point for modelling the structure of the network. Hereby in order to better understand the inside of such a stochastic system, we need to find a mathematical model which can describe its behaviour and to estimate its model parameters. There are different mathematical models which enable us to represent the complex biological systems by linear and nonlinear ways [3–6] and to infer the model parameters stochastically [7]. Among modelling approaches, the nonlinear ones, which can be listed as the Langevin [8], diffusion approximation [9, 10], and inhomogeneous poisson process [10], can capture the actual randomness in the system by using the chemical master equation [11]. In these approaches, the first two models explain the system by the same expression in the sense that the former converts a deterministic model to stochastic one by adding a noise term in the equation [8]. Then it solves this extended model, that is, noise-added deterministic model, by Itô or Stratonovich integrals [11, 12]. On the other hand, the latter converts a stochastic expression to the differential equations via the Fokker-Planck equation [11]. The final approach, which is the inhomogeneous poisson process model, is based on the Gillespie technique which can exactly simulate the biological network stochastically [13]. However, since it is computationally demanding in inference of the model parameters, currently it is implemented for toy systems [10]. Therefore in this study we consider to implement the nonlinear diffusion approximation to model a realistically large MAPK/ERK pathway as it enables us to estimate the parameters via convex optimization techniques. In particular, we perform the Markov chain Monte Carlo methods for the inference of the NP hard problem.Indeed for the purpose of estimation, there are several alternative methods. One of them is the application of the parallel computing [14]. In a Bayesian setting the parallelization technique can be effective in simple situations when the model parameters are conditionally independent on one another. Whereas it is not useful when the blocks of dependent variables are large [10, 14], like in our biochemical system. Another recent alternative is the method discussed in [6]. They propose two approaches. The first approach is to use the reversible jump methods along with Metropolis-Hastings acceptance schemes and the second is to implement the block updating methods based on the use of Poisson process approximations and random walk schemes. From their analysis in the Lotka-Volterra model [10], they show that both methods are promising in efficiency; however, further extensions for the implementations in complex networks are needed. For the inference of the complex network dynamics, [15] suggests the application of the deterministic Michaelis-Menten model whose measurements possess independent log-normally distributed location and scale parameters. The estimates of model parameters are computed by the discretized maximum likelihood method via a piecewise constant step function in the small time interval. On the other hand, [16] uses a regression model whose regression coefficients are different for all genes for the given transcription factor activities (TFAs). The inference of the parameters is conducted by optimizing the likelihood of observations via scaled conjugate gradient algorithm. In this estimation as the number of model parameters becomes very large due to the assumption of different TFAs for each gene, a dimension reduction technique is implemented by putting a Gaussian prior on the transcription factors as the constraints in the marginalization of the likelihood function. Apart from these techniques there are some other approaches which are based mainly on the ordinary differential equations (ODEs) and hybrid methods which solve the nonlinear structure of the equilibrium state in ODE by adding noise terms [8, 17]. For instance, [18] implements ODE methods in conjunction with a Bayesian setting for the inference of a simple MAPK system. Reference [19] estimates the model parameters of the JAK2-STAT5 pathway by ODEs whose model is represented by power-law terms, resulting in noninteger kinetic orders. In that study these kinetic orders, which are the model parameters in this case, are solved via the genetic algorithm. The aim of this paper is to check the global structure of the MAPK/ERK system on western blot time-course data of a subset of proteins [20] by stochastic approach. As stated beforehand, we consider a diffusion approximation technique which enables to deal with the realistic complexity while maintaining the computational efficiency. In fact, the diffusion method uses continuous time and requires infinite information to describe a finite sample path, when the measurements are discrete and are only partially observed. Thus, we implement a discretized version of the diffusion approximation called the Euler-Maruyama approximation. This discretizes time, but keeps the continuous structure of measurements. The details of this method are described in [21]. But in this paper, different from the study in [21], we consider that the randomness in the model is also caused by the measurement error in the western-blotting dataset. Therefore as described in Section 3, we construct a more realistic model for the inference and change the updating stages of MCMC computation accordingly. Moreover in the calculation we use very limited observations for both the time points and species which enables us to check both the applicability of our proposal mathematical model and inference algorithm. Furthermore, in this study we suggest biologically plausible descriptions for the pathway and select the best-fitted model for our dataset by model-selection criteria and in the end we evaluate whether the selected model can really validate the biological findings about the system. Hereby in Section 2 we describe alternative models for the pathway structure. Section 3 outlines the formulation of the Euler-Maruyama approximation by means of the hierarchical model and gives details about the MCMC updates of the system. In Section 4.1 we describe the available western blot dataset and in Section 4.2 we present the results and compare the models via several information criteria. Finally, in Section 5, we conclude the outputs.


## 2. System and Methods

Due to its importance, the MAPK pathway has been studied intensively. Hence there are a number of sources that provide knowledge about this system. Most of this information is qualitative, and only a small part of it has been captured in an explicit set of reactions. Here we combine this qualitative information to present the steps of the activation of the pathway as a list of quasi reactions. A novelty of our approach is the usage of multiple state variables to deal with molecules for which different localizations in the cell are an intricate part of the dynamic process and to handle proteins that have different binding sites and various phosphorylation states. For expressing the translocation of proteins to the membrane, we use the notationm. On the other hand, the different levels of the phosphorylation are denoted by distinct abbreviations in the sense that p or p1 show the monophosphorylation, whereas p2 indicates the double-phosphorylation of a protein. For instance, in the following set of equations:(a)
Raf+PP2A→c1Raf.I+PP2A;(b)
Raf.I+Ras.GTP→c2Raf.Im+Ras.GTP;(c)
MEK.p2+ERK→c3MEK.p2+ERK.p1;(d)
MEK.p2+ERK.p1→c4MEK.p2+ERK.p2;we describe the interaction between PP2A and inactive Raf proteins by Reaction (a). This expression implies that the PP2A protein takes away the inhibitory phosphorylate of the inactive Raf (Raf) which enables Raf to activate. The underlying inactive and monophosphorylated Raf (Raf.I) is translocated by the recruitment of the active Ras (Ras.GTP) from the cytosol to the cell membrane (Raf.Im) as denoted by Reaction (b). The activation of the ERK protein is shown by Reactions (c) and (d). Here the inactive ERK protein is double-phosphorylated via the active MEK (MEK.p2) with two steps, where ERK.p1 and ERK.p2 illustrate the mono- and double-phosphorylation of ERK proteins in the corresponding reactions. Using this representation we describe the pathway via 51 species involving 34 major proteins. Similarly, we use 94 reactions where 65 of them reflect changes in activities and translocations of species and the rest shows their degradations after dissociation. In biochemical reactions the probability per unit time of each event is stated by a stochastic rate constant. Here those constants are denoted as ci’s (i=1,…,94) which are the parameters to be estimated.For inferringci, we implement MCMC algorithms within a Bayesian framework. These algorithms enable us to deal with the challenges of genomic measurements, that is, missing data and a limited knowledge of the complete network structure. Moreover, they are able to include the uncertainty coming from various sources of variations [18].We define various alternative models for the MAPK system by either combining several reactions in a single equation, hereby simplifying the system via reducing the number of total reactions, or by separating the pathway into small parts in such a way that the resulting modules will still be able to explain the biological relationships in the experiments. We suggest 5 models by using the listed totals below. But when generating these totals, we accept several assumptions by considering that the data are gathered by western blot: (i) the weight of protein complexes is heavier than the weight of single proteins, whereas (ii) single proteins have the same weight regardless whether they are mono/double-phosphorylated:(a)
TotalRas=Ras.GDP+Ras.GTP,(b)
TotalRaf=Raf+Raf.I+Raf.Im+Raf.Am,(c)
TotalMEK=MEK+MEKS+MEKF+MEK.p2,(d)
TotalERK=ERK+ERK.p2for Models 1 and 2,(e)
TotalERK=ERK+ERK.p1+ERK.p2for Models 3, 4, and 5,where Ras.GDP and Ras.GTP represent inactive and active Ras substrates, respectively.MEKS and MEKF denote MEK proteins phosphorylated by PAK and ERK proteins, in the given order. Raf.I and Raf.Im designate the inactive Raf in the cytosol and near the cell membrane, respectively, and finally Raf.Am indicates the active Raf near the cell membrane. Due to the available measured totals we face a convolution problem, for example, with respect to Raf, Raf.I, Raf.Im, MEK, MEKS, and MEKF proteins in all models. Models 3, 4, and 5 also have this problem with ERK and ERK.p1 proteins. To use the information from those totals, we need additional constraints in inference. We assume that for observed time points the number of molecules of MEK proteins should be equal to MEK=TotalMEK-(MEKS+MEKF+MEK.p2), where MEKS and MEKF are both positive. Also for most Raf measurements we observe negative values after the substraction. Therefore, we completely omit the information on the Total Raf from all models. All suggested models are listed as below and their relationships are simply drawn in Figure 1. On the other hand, the complete list of reactions for each model with the description of their associated proteins can be found in Appendices B–G.Figure 1
Simple representation of the relationship between nested structure of the five suggested MAPK/ERK pathway models.Model 1.
According to the current biological theory of the MAPK pathway [20], it is suggested that RKIP and ERK proteins are the two main components which regulate the inhibitory control of the system. We assume that this control can be directed without RKIP since the complete activity, including the inhibitory control, of the pathway is still possible without RKIP terms as long as ERK and its associated complexes exist in the system. This assumption enables the reduction in the number of substrates from 51 to 41 without altering the main activation of the system.
Moreover, seeing that in biochemical reactions, the protein degradations are much slower than the time periods during which biochemical activation and deactivation processes take place, we consider that protein degradations are insignificant and can be canceled from the reaction list. We also assume that some reactions can be combined into a single reaction if they are executed almost simultaneously or if they introduce another variable which would have no other function than to react with proteins in the given reaction. For instance, the equationsMEK.p2+ERK→MEK.p2+ERK.p1    andMEK.p2+ERK.p1→MEK.p2+ERK.p2 are summarized as MEK.p2+ERK→MEK.p2+ERK.p2 resulting in the cancelation of ERK.p1 proteins. By applying this simplification, we describe the system via 39 rather than 47 reactions, and 38 rather than 41 substrates in which 20 of them are linearly independent. In our estimation, we particularly eliminate merely nonobserved substrates because of their linear dependencies on other species (see Section 3), thereby all time-course observations are used in calculations. So the estimation is conducted by 6 observed (Ras.GDP, Ras.GTP, Raf.Am, MEK.p2, ERK, and ERK.p2) and 14 unobserved substrates.Model 2.
It is the extension of Model1 by including the degradation of the EGF receptor (EGFR). The EGFR degradation is the direct result of the MAPK activation by the internalization into vesicles of this receptor, thereby can be important for the steady-state description of the system. Thus, the inference is based on 40 reactions and 38 substrates in which 21 of them are linearly independent and the same 6 substrates are used as observed.Model 3.
Model1 is extended by describing all reactions rather than summarizing reaction groups into a single equation. Apart from those changes, the model still excludes all degradations and all reactions which have RKIP proteins or their complexes. In Model 3 the pathway is represented by 47 reactions and 41 substrates where 26 of them are linearly independent and 5 (Ras.GDP, Ras.GTP, Raf.Am, MEK.p2, and ERK.p2) among 26 substrates have observed measurements. All observed substrates are included in computations and the elimination of substrates due to linear dependence is conducted among the nonobserved ones. But different from previous models, the explicit measurements of ERK proteins are lost because of the convolution between ERK and ERK.p1 proteins.Model 4.
It is an extension of Model3 in the sense that the EGFR degradation is included into the system because of the effect of the EGFR internalization as taken in Model 2. This results in 48 reactions and 41 substrates. In the calculation we use 27 linearly independent substrates where the common 5 substrates have real time-course measurements.Model 5.
This is the largest model such that the pathway includes all types of RKIP proteins, accordingly its inhibitory control is regulated via both ERK and RKIP species. Moreover it is assumed that the EGFR degradation is essential for the best MAPK description and the system does not use any kind of summary reactions. Hence, the inference is conducted by 66 reactions and 51 substrates. But due to the elimination of linearly dependent ones, the computation is done via 29 unobserved and 5 common observed proteins.


## 3. Algorithm

Under the assumption that the probability distribution of the number of moleculesY in each species at time t, P(Y,t), is continuous, the stochastic model can be converted to a differential equation model. Indeed, P(Y,t) can be further expanded via a Taylor series expansion and the change in states of each species at t is found by a Fokker-Planck approach, in which a correlated noise term describes the stochastic behaviour of the model over and above a deterministic drift term [11]. Under this condition a nonlinear Fokker-Planck equation, which is converted from the continuous approximation of the chemical master equation, can be written in terms of infinitesimal mean and second-moment rates of the jump in states. By applying the Itô diffusion to the underlying Fokker-Planck equation, we get a diffusion formulation of the system. In a diffusion process, a finitely observed sample path is strictly speaking intractable. But under the assumption of discrete jumps at a large set of discrete time points, we can employ the Euler-Maruyama approximation of the diffusion process [10, 22, 23], which is one of the common techniques in stochastic modelling where the optimal values of the model parameters are found from the convergence distribution of the MCMC runs. Hereby the underlying model is described as follows:
(3.1)ΔYt=μ(Yt,Θ)Δt+β1/2(Yt,Θ)ΔWt,
in which μ(Y,Θ)=V′h(Y,Θ) and β(Y,Θ)=V′diag{h(Y,Θ)}V are mean, or drift, and variance, or diffusion, matrices, respectively, both depending on the state Y=(Y1,…,Yn) at time t and the parameter vector Θ=(c1,…,cr)′ explicitly. n is the total number of substrates and r states the total number of reactions in the system. ΔWt is an n-dimensional independent identically distributed Brownian random vector over time Δt. V is the net effect matrix which is the difference between the stoichiometric coefficients of products and reactants of each reaction. Finally, h(Y,Θ) indicates the hazard, also called the rate law, of the reaction [10]. In these expressions, the notation (′) denotes the transpose of the underlying matrix or vector.In the estimation of the reaction rates of the MAPK system, we consider that the state matrixY is composed of both observed species X and unobserved or missing species Z at given time point t, hereby Yt≡(Xt,Zt)′ where each Y at t consists of totally n species as previously stated. To sample the missing states given the parameters and the parameters given all states, we use Metropolis-within-Gibbs steps. This algorithm is suitable for the cases in which the Gibbs and Metropolis-Hasting algorithms can be used in a cycle [24]. In the MCMC, furthermore, we apply data augmentation for the nonobserved states in the Euler process by putting latent states within each pair of time point. In this way, the strong bias in estimators caused by the large time steps in the Euler approximation is decreased in exchange for some additional computational complexity. Hereby we augment the state matrix Y, whose observations are not evenly spaced in time, in such a way as to make the resulting full state equally spaced in time. We consider two different regimes with more and less dense states. Accordingly, we add states in every 5 minutes and 10 minutes between each pair of totally 8 observed time points (Table 1) resulting in 37 and 77 time points, respectively, in observation matrices of all suggested models.Table 1
Western blotting time-course dataset for the MAPK/ERK pathway [20].
Time (min)
Ras.GTP
Total Ras
Raf.Am
Total Raf
0
10
175
84
231
20
67
189
307
239
40
108
278
282
243
60
60
328
441
249
120
202
252
227
221
180
157
374
321
311
240
148
348
244
266
360
156
295
289
321
Time (min)
MEK.p2
Total MEK
ERK.p2
Total ERK
0
19
217
28
694
20
186
202
133
682
40
205
210
191
682
60
152
268
151
582
120
157
250
82
673
180
137
238
167
710
240
141
289
83
572
360
126
240
115
460Finally, in our model we assume that each observation at timet possesses a measurement error, thereby, the model uncertainty originates from not only the stochasticity of the protein interactions, but also the variabilities from the measurement error. The application of this assumption in financial volatility models can be found in [25]. On the other hand, [5] use this assumption for the inference of a simple model of prokaryotic autoregulation. Here this assumption is inserted in our complex MAPK model by defining every observation at t under measurement error, Wt, via Wt=Xt+σm, where σm denotes the measurement error coming from normal distribution with mean 0 and variance, also called the tuning parameter, φ2 for all time and observed substrates, that is, σm~N(0,φ2). In the sampling, since the tuning parameter φ2 has an effect on the acceptance rate of σm, φ2 is adjusted adaptively during the first 10,000 iterations of the burn-in period in the MCMC algorithm. More details about this sampling and the augmentation of Y can be found in [21].In inference of the parametersΘ=(c1,…,cr,σm), the positivity of rate constants ci (i=1,…,r) is maintained by independent exponential priors with rate 1. For the prior of the measurement error σm, we select the inverse gamma with scale and location parameters 1 due to its conjugality of the normal distribution. Both densities produce flat and heavily tailed prior distributions for the elements of Θ. In the update of the system, we use block updates to sample the reaction rates ci. In our block scheme we divide the r-dimensional vector of c into small and equally sized groups with a dimension d and then simulate each d-dimensional group sequentially. The size of the block update dimension d in each model is varied. We set d to 3, 4, and 5 for Models 1 and 2, Models 3 and 4, and Model 5, respectively. After the update of rate constants, we update σm.When we renewci (i=1,…,r), we begin to update Y column by column whose entries are composed of partially observed and completely augmented data. The candidate values of each Y column are generated from distinct transition kernels via the Metropolis-within-Gibbs algorithm. We include the measurement error in the update of the partially observed states Wt from a normal distribution with mean Xt and variance σm via Wt~N(Xt,σm)·Xt denotes the observed values of the substrates at the given time t as described in advance. On the other side, the proposal of the completely augmented states is produced as stated in [9]. In Appendix A we present the derivations of candidate generators for partially observed states as an example of transition probabilities. At the end of the update Y, the convergence of the chain is controlled and the algorithm is repeated until the convergence is attained.During the application of these MCMC methods to the five MAPK pathway models, we suggest an updating scheme which is able to deal with any kind of dependency coming from the complexity of our system in the sampling process. We mainly divide the cause of dependency into two parts, namelystructural and incidental dependence. The former is originated from the linear dependence of the substrates in the net effect matrix V′V. We eliminate these problematic species when starting the algorithm so that there is no singularity problem in the calculation of the likelihood. The latter, on the other hand, is caused by the singularity of β(Y,Θ)=V′diag{h(Y,Θ)}V or any numeric problems resulting from h(Y,Θ) where h(Y,Θ) is the n-dimensional hazard vector of the system for a given Y and Θ, and diag{h(Y,Θ)} represents the n×n diagonal matrix of h(Y,Θ). To work with the incidental dependency we check whether a candidate value for either missing states of Y or Θ causes a singularity in the diffusion if it were accepted. If it preserves the nonsingularity, the acceptance probability is computed and the system is updated in the usual way. Otherwise, the candidate value is rejected even before calculating the acceptance probability. The details about this MCMC scheme can be found in [21].


## 4. Implementation

### 4.1. Description of the Experimental Data

We use the experimental measurements of the protein levels and activities of Ras, Raf, MEK, and ERK proteins from mitogen stimulated COS1 green monkey kidney cells. These cells can easily be propagated in the cell culture and are commonly used for biochemical experiments. The cells and their crude protein extracts contain many ten thousands of different proteins such that the proteins of interest constitute only a small proportion of the mixture. Hence the cell extracts usually need to be separated in order to be able to visualize the protein of interest. In collecting the experimental protein data from the underlying cell extracts, the technique of thewestern blotting, the immunoprecipitation, and the pulldown are widely utilized. If the proteins in a sample of the dissolved cell are detected by separating proteins according to their weights, the method is called the western blotting [26]. This method is used here to measure MEK and ERK activities. If a further purification is required, in order to measure enzymatic activities, the immunoprecipitation technique is used. So the protein of interest from the crude is extracted and its enzymatic activities are assayed by the incubation with relevant substrates in vitro. This is the way how the Raf activation is measured. A pulldown assay is a variation of the immunoprecipitation in the sense that the purification of pulldown is done by means of an assay specific to the binding side rather than a substance. In our data the measurements of Ras proteins are collected by this technique [27].In our dataset initially the expression levels of these proteins which do not change over the time of the mitogen stimulation are examined as the control expressions and then the measured activities are reflected as the relative changes between untreated and mitogen-stimulated cells. The ODE method is successful in working with such changes in concentrations but they are insufficient to present the stochastic manner of biological activations [28, 29]. On the contrary, the stochastic methods enable to take these effects into account by working with the number of molecules. By assuming that the measurements are proportional to the true underlying number of protein molecules, we arbitrarily set the lowest intensity in the order to 10 molecules and extrapolate all other molecules from that proportion as shown in Table 1. On the other hand, the actual time-course data without extrapolation are represented in Table 5 and are plotted in Figure 5 in Appendix.


### 4.2. Application to the Experimental Data

In inference of all models, the total number of iterations is chosen as 100,000 in which the first 85,000 are taken as the burn-in. The full list of reactions and concerning estimations are presented in AppendicesB–F. From these results, first we check how strong the prior distributions of parameters affect their probability distributions. For the practical purpose we choose the Kolmogorov-Smirnov test in calculations, thereby compare MCMC outputs of rate constants and measurement error to random variables generated from exponential distribution with rate 1 and inverse gamma with parameters 1, respectively. The test statistics indicate significant differences for all models, empirically implying that the selected prior densities in fact cause no effect on the corresponding posterior densities. Then we compare the estimates found at Δt=5 and Δt=10. Using PAM clustering [30] with k=5 classes (referring to very slow, slow, moderately slow, moderately fast, very fast) on the parameter estimates, we evaluate their stabilities. Table 2 shows that the estimated values are very stable between the two approximations. Moreover, from the analysis we see that although the majority of model parameters possess good mixing in the sense that the acceptance ratios are not less than our lower bound of 0.05, some of them own low mixing property such as 0.02 or 0.03. We consider that the sparsity of our dataset, high dependencies between species, and the existence of latent states cause less-likely candidate values in MCMC runs although the adaptive sampling plan during the first 10,000 MCMC runs enables us to get good tuning parameters. Furthermore, from the comparison of the computational demands under both Δt’s, it is observed that the increase of the augmented states by going from Δt=10 to Δt=5 significantly raises the computational time (Table 3). On the other hand, from the comparison of estimates, it is found that the correlations between MCMC runs are slightly higher in Δt=5 than in Δt=10. Moreover we get better mixing in some of the estimates under Δt=5 which is consistent with [9] since the augmentation in latent states leads to higher dependencies between states. But we do not observe a complete improvement for all rates under each model. We consider that the reason is relevant to the structure of our observed time points. Because for our western blotting data, the number of augmented states rises after the observed time point t=60min (12 latent states from t=60 to t=180min and 24 latent states from t=240 to t=360min), hereby, the dependency between the parameters and those states increases. On the other side, since the number of augmented states is not large from t=0 to t=60min, the underlying augmentation does not very much results in worsening the dependency among the updates of these particular states and estimates (4 latent states from t=0 to t=60min). Finally, it is seen that the estimated measurement errors are very big with respect to the observations in Table 1. We consider that similar to the estimates of reaction rates, the high correlation between the species as well as the sparsity of the observed time points and observed species lead to low mixing property, thereby, high measurement errors.Table 2
Stability of the PAM clustering of alternative MAPK/ERK pathway models via the posterior means of model parameters forΔt=5 and Δt=10 with 5 clusters. The estimates are based on 15,000 MCMC runs after 85,000 burn-in runs.
Model
Number of parameters
Number of changes fromΔt=5 to Δt=10
Model1
39
1
Model2
40
1
Model3
47
5
Model4
48
0
Model5
66
1Table 3
Computational time (where h, min, and sec stand for hour, minute, and second, resp.) for estimating model parameters via 100,000 MCMC runs of all models whenΔt=5 and Δt=10 in R on the 3.00 GHz Dual Core Xeon processor with a single-trade application.
Model
Δ
t
=
5
Δ
t
=
10
Real time
Real time
1
42 h 35 min 36 sec
21 h
2
45 h 23 min 54 sec
22 h 53 min 40 sec
3
52 h 26 min 43 sec
27 h 35 min 17 sec
4
56 h 01 min 41 sec
28 h 52 min 05 sec
5
91 h 46 min 42 sec
26 h 57 min 21 secAs a result of all these outcomes we think that though the current data are not very reliable for finding the possible reaction rates of the MAPK system, the estimates can be still used for the procedure of the model selection since the inference of all models is conducted under the same challenges and the main interest from the estimation results is to investigate whether any proposal model for the MAPK pathway (Section2) suggests a better fit for the available western blotting dataset (Table 1). Although there are a number of model selection methods in literature, we choose the Bayesian information criteria (BIC) and deviance information criteria (DIC) as the model selection methods since they give a distance of the data to each of the alternative models which are not generalized from others, and enable us to compare all models even if all of which are incorrect in a biological sense [31].In inference, the number of associated substrates having the structural dependence is different for each model. But considering our sparse western blotting measurements, we particularly order the species for every model in such a way that none of the observed substrate presented in Table1 is discarded due to the structural dependency. In this way, we can execute the complete measurements in Table 1 in our estimation, thereby implement the comparison of both BIC and DIC via all data given in Table 1 in each alternative model. In our analysis although we observe slightly high autocorrelation under Δt=5 than the estimates under Δt=10, we still choose the results from Δt=5 for the model selection seeing that the estimates of the associated measurement errors are lower that the ones under Δt=10. Table 4 summarizes the test results based on MCMC outputs sampled from the convergent parts of the chains which are taken as the last 15,000 runs for each model. In order to select samples for the comparison, we check the autocorrelation functions (ACF) of each substrate per model and take 500 thinned samples whose MCMC outputs have relatively less correlations within each other. From the tabulated BIC and DIC values it is clear that Model 2 is the best fitting model for our western blot data. This model suggests that the system does not include any kinds of RKIP proteins and its components, and can be described by summarizing simultaneous reactions into a single equation. Moreover, it indicates that the degradation of EGFR is essential in the description of this system.Table 4
Results of BIC and DIC criteria for the alternative MAPK/ERK pathway models based on 500 MCMC samples whenΔt=5 by data in Table 1.
Model1
Model2
Model3
Model4
Model5
BIC
183622.8
28450.14
76446.05
41631.55
136870.2
DIC
325701
46548.61
128999.5
66110.5
186086.0Table 5
Western blotting time-course dataset for the MAPK pathway collected from the School for Cancer Studies in Glasgow University.
Time (min)
Ras.GTP
Total Ras
Raf.Am
Total Raf
0
35618
623520
300095
824192
20
240222
672847
1092111
850161
40
385132
989878
1005269
865119
60
212275
1169001
1570648
887820
120
718095
898817
810234
785545
180
559659
1332240
1142328
1105969
240
525533
1240655
867336
948456
360
555147
1049890
1028401
1141482
Time (min)
MEK.p2
Total MEK
ERK.p2
Total ERK
0
67914
772901
98683
2470824
20
662655
718812
472934
2430305
40
729792
748227
678992
2428199
60
539726
954104
536933
2071863
120
559206
889207
290295
2396161
180
487896
846376
594756
2529436
240
501463
1030196
295292
2038472
360
449943
856206
408216
1636536In a biological sense Model2 implies that the inhibited control of the Raf-MEK pathway is directed via ERK, rather than ERK and RKIP, proteins by phosphorylating MEK at T292 binding site or by phosphorylating SOS, as long as the EGR receptor triggers the activation of the pathway. Moreover, the chains of reactions which produce, particularly, Shc-Grb2-SOSm, Grb2-SOS, ERK.p2, and Ras.GTP are very quick in the sense that Raf.A-Ras.GTPm, ERK.p1, and Shcm proteins, which are used in the production of those proteins, cannot stay longer in the cell to count them as separate species. On the other hand although the knowledge about the MAPK system claims the existence of all kinds of RKIP proteins in the system [32–34], Model 2 describes the system without any RKIP species.Once the selection of the best-fitted model we compare it with the proposal models of the system in [35] which are suggested by [36–38]. The main differences among these models are the function of the EGF receptor during the activation, the existence of Shc substrates and negative-feedback loops. Among these alternative models, Model 2 supports the model of [36]. The main feature of that model is the existence of the Ras activation with Shc and without Shc species, the degradation of EGRF, and core cascade of ERK via Raf and MEK. Moreover, different from this model, Model 2 includes the negative feedback phosphorylation of SOS by active ERK and defines the system with less species and reactions. Schoeberl et al. [36] model describes the pathway with 94 reactions and 125 reactions. On the other hand to the best of our knowledge since none of this model is validated by a real time-course dataset within a stochastic modelling, there have not yet any wet-data/approaches which can verify Model 2 or other suggested models.On the other hand, if we only evaluate the estimated values of reaction rates of Model2, we see that the dissociation type of reactions is the fast reactions and among other dissociations, the ones within the inactive Raf and active Ras complex (reaction 15 in Model 2) as well as within the active ERK, RSK, and transcription factor (reaction 39 in Model 2) have the highest speed. Whereas the reactions which regulate the Raf-MEK cascade by the inhibition of ERK (reactions 20 and 23 in Model 2) occur slowly. Figure 2 is shown as example from the posterior distributions of some reaction rates and the trace plots of their MCMC outputs. Figure 3, on the one hand, briefly describes Model 2. In this figure the thickness of the interaction indicates the speed of the associated reactions.Examples of the probability distribution (a, b, and c) and trace plots (d, e, and f) of reaction rates 2 (which indicates the simultaneous recruitment of Shc and Grb2-SOS complex from the cytosol to the cell membrane by the recruitment of EGFR (a and d)), 16 (which shows the activation of MEK proteins by active Raf (b and e)), and 38 (which refers to the dissociation of active ERK and RSK complex (c and f)) of Model2 under Δt=5 after 85,000 burn-in runs.
(a)
(b)
(c)
(d)
(e)
(f)Figure 3
Simple representation of the MAPK/ERK pathway via Model2 under Δt=5. The thickness of the arrow indicates the speed of the interaction. The thick arrows represent fast interactions (estimated reaction rates have power 10-1), only simple arrows describe moderately slow reactions (estimated reaction rates have power 10-2 or 10-3), and finally dotted lines show very slow reactions among substrates (estimated reaction rates have power 10-4 or 10-5).Figure 4
Simple representation of the structure of the MAPK/ERK pathway taken from Figure 2 in the study of [21].Figure 5
Plot of the western blotting data of the MAPK pathway observed at 8 time points according to the given molecular populations in Table5.On the other hand, [39] shows that the single element update of states causes poor mixing due to the high correlation within the latent data. Thereby as an extension of this study, we can investigate the possibility of block updates of missing states [40, 41] to get better mixing. We consider that this suggested approach improves the current approximation technique in the class of stochastic methods where the optimal solution of the system cannot be found globally. Moreover, as shown in Figure 2, the curves may also imply fractal time series [42, 43] in which the system can be also modelled via the fractional Brownian motion (FBM). In FBM, the density of the time series has heavy-tailed distribution and indicates slowly decreasing autocorrelation functions with power spectrum of power-law type that can be observed as different cell signalling processes [44] as well as earthquake modelling and financial time series analysis. The power-law spectrum is seen if the process shows the self-similarity, that is, a time segment in the process presents the same behaviour as any segments of other time scales, and stationarity of its increments in diffusion model [45, 46]. If the self-similarity is found locally, rather than globally as in FBM, the system can be modelled by a more general modelling approach, called the multifractional Brownian motion [45]. On the other side if we deal with the change in increment of the Gaussian process, the system can be described by the fractional Gaussian model (FGM) [47–50]. Hereby as the extension of the diffusion modelling, the high correlation in the system can be solved by both FBM and FGM models that are specifically designed for such short- or long-correlated time series datasets.


## 4.1. Description of the Experimental Data

We use the experimental measurements of the protein levels and activities of Ras, Raf, MEK, and ERK proteins from mitogen stimulated COS1 green monkey kidney cells. These cells can easily be propagated in the cell culture and are commonly used for biochemical experiments. The cells and their crude protein extracts contain many ten thousands of different proteins such that the proteins of interest constitute only a small proportion of the mixture. Hence the cell extracts usually need to be separated in order to be able to visualize the protein of interest. In collecting the experimental protein data from the underlying cell extracts, the technique of thewestern blotting, the immunoprecipitation, and the pulldown are widely utilized. If the proteins in a sample of the dissolved cell are detected by separating proteins according to their weights, the method is called the western blotting [26]. This method is used here to measure MEK and ERK activities. If a further purification is required, in order to measure enzymatic activities, the immunoprecipitation technique is used. So the protein of interest from the crude is extracted and its enzymatic activities are assayed by the incubation with relevant substrates in vitro. This is the way how the Raf activation is measured. A pulldown assay is a variation of the immunoprecipitation in the sense that the purification of pulldown is done by means of an assay specific to the binding side rather than a substance. In our data the measurements of Ras proteins are collected by this technique [27].In our dataset initially the expression levels of these proteins which do not change over the time of the mitogen stimulation are examined as the control expressions and then the measured activities are reflected as the relative changes between untreated and mitogen-stimulated cells. The ODE method is successful in working with such changes in concentrations but they are insufficient to present the stochastic manner of biological activations [28, 29]. On the contrary, the stochastic methods enable to take these effects into account by working with the number of molecules. By assuming that the measurements are proportional to the true underlying number of protein molecules, we arbitrarily set the lowest intensity in the order to 10 molecules and extrapolate all other molecules from that proportion as shown in Table 1. On the other hand, the actual time-course data without extrapolation are represented in Table 5 and are plotted in Figure 5 in Appendix.


## 4.2. Application to the Experimental Data

In inference of all models, the total number of iterations is chosen as 100,000 in which the first 85,000 are taken as the burn-in. The full list of reactions and concerning estimations are presented in AppendicesB–F. From these results, first we check how strong the prior distributions of parameters affect their probability distributions. For the practical purpose we choose the Kolmogorov-Smirnov test in calculations, thereby compare MCMC outputs of rate constants and measurement error to random variables generated from exponential distribution with rate 1 and inverse gamma with parameters 1, respectively. The test statistics indicate significant differences for all models, empirically implying that the selected prior densities in fact cause no effect on the corresponding posterior densities. Then we compare the estimates found at Δt=5 and Δt=10. Using PAM clustering [30] with k=5 classes (referring to very slow, slow, moderately slow, moderately fast, very fast) on the parameter estimates, we evaluate their stabilities. Table 2 shows that the estimated values are very stable between the two approximations. Moreover, from the analysis we see that although the majority of model parameters possess good mixing in the sense that the acceptance ratios are not less than our lower bound of 0.05, some of them own low mixing property such as 0.02 or 0.03. We consider that the sparsity of our dataset, high dependencies between species, and the existence of latent states cause less-likely candidate values in MCMC runs although the adaptive sampling plan during the first 10,000 MCMC runs enables us to get good tuning parameters. Furthermore, from the comparison of the computational demands under both Δt’s, it is observed that the increase of the augmented states by going from Δt=10 to Δt=5 significantly raises the computational time (Table 3). On the other hand, from the comparison of estimates, it is found that the correlations between MCMC runs are slightly higher in Δt=5 than in Δt=10. Moreover we get better mixing in some of the estimates under Δt=5 which is consistent with [9] since the augmentation in latent states leads to higher dependencies between states. But we do not observe a complete improvement for all rates under each model. We consider that the reason is relevant to the structure of our observed time points. Because for our western blotting data, the number of augmented states rises after the observed time point t=60min (12 latent states from t=60 to t=180min and 24 latent states from t=240 to t=360min), hereby, the dependency between the parameters and those states increases. On the other side, since the number of augmented states is not large from t=0 to t=60min, the underlying augmentation does not very much results in worsening the dependency among the updates of these particular states and estimates (4 latent states from t=0 to t=60min). Finally, it is seen that the estimated measurement errors are very big with respect to the observations in Table 1. We consider that similar to the estimates of reaction rates, the high correlation between the species as well as the sparsity of the observed time points and observed species lead to low mixing property, thereby, high measurement errors.Table 2
Stability of the PAM clustering of alternative MAPK/ERK pathway models via the posterior means of model parameters forΔt=5 and Δt=10 with 5 clusters. The estimates are based on 15,000 MCMC runs after 85,000 burn-in runs.
Model
Number of parameters
Number of changes fromΔt=5 to Δt=10
Model1
39
1
Model2
40
1
Model3
47
5
Model4
48
0
Model5
66
1Table 3
Computational time (where h, min, and sec stand for hour, minute, and second, resp.) for estimating model parameters via 100,000 MCMC runs of all models whenΔt=5 and Δt=10 in R on the 3.00 GHz Dual Core Xeon processor with a single-trade application.
Model
Δ
t
=
5
Δ
t
=
10
Real time
Real time
1
42 h 35 min 36 sec
21 h
2
45 h 23 min 54 sec
22 h 53 min 40 sec
3
52 h 26 min 43 sec
27 h 35 min 17 sec
4
56 h 01 min 41 sec
28 h 52 min 05 sec
5
91 h 46 min 42 sec
26 h 57 min 21 secAs a result of all these outcomes we think that though the current data are not very reliable for finding the possible reaction rates of the MAPK system, the estimates can be still used for the procedure of the model selection since the inference of all models is conducted under the same challenges and the main interest from the estimation results is to investigate whether any proposal model for the MAPK pathway (Section2) suggests a better fit for the available western blotting dataset (Table 1). Although there are a number of model selection methods in literature, we choose the Bayesian information criteria (BIC) and deviance information criteria (DIC) as the model selection methods since they give a distance of the data to each of the alternative models which are not generalized from others, and enable us to compare all models even if all of which are incorrect in a biological sense [31].In inference, the number of associated substrates having the structural dependence is different for each model. But considering our sparse western blotting measurements, we particularly order the species for every model in such a way that none of the observed substrate presented in Table1 is discarded due to the structural dependency. In this way, we can execute the complete measurements in Table 1 in our estimation, thereby implement the comparison of both BIC and DIC via all data given in Table 1 in each alternative model. In our analysis although we observe slightly high autocorrelation under Δt=5 than the estimates under Δt=10, we still choose the results from Δt=5 for the model selection seeing that the estimates of the associated measurement errors are lower that the ones under Δt=10. Table 4 summarizes the test results based on MCMC outputs sampled from the convergent parts of the chains which are taken as the last 15,000 runs for each model. In order to select samples for the comparison, we check the autocorrelation functions (ACF) of each substrate per model and take 500 thinned samples whose MCMC outputs have relatively less correlations within each other. From the tabulated BIC and DIC values it is clear that Model 2 is the best fitting model for our western blot data. This model suggests that the system does not include any kinds of RKIP proteins and its components, and can be described by summarizing simultaneous reactions into a single equation. Moreover, it indicates that the degradation of EGFR is essential in the description of this system.Table 4
Results of BIC and DIC criteria for the alternative MAPK/ERK pathway models based on 500 MCMC samples whenΔt=5 by data in Table 1.
Model1
Model2
Model3
Model4
Model5
BIC
183622.8
28450.14
76446.05
41631.55
136870.2
DIC
325701
46548.61
128999.5
66110.5
186086.0Table 5
Western blotting time-course dataset for the MAPK pathway collected from the School for Cancer Studies in Glasgow University.
Time (min)
Ras.GTP
Total Ras
Raf.Am
Total Raf
0
35618
623520
300095
824192
20
240222
672847
1092111
850161
40
385132
989878
1005269
865119
60
212275
1169001
1570648
887820
120
718095
898817
810234
785545
180
559659
1332240
1142328
1105969
240
525533
1240655
867336
948456
360
555147
1049890
1028401
1141482
Time (min)
MEK.p2
Total MEK
ERK.p2
Total ERK
0
67914
772901
98683
2470824
20
662655
718812
472934
2430305
40
729792
748227
678992
2428199
60
539726
954104
536933
2071863
120
559206
889207
290295
2396161
180
487896
846376
594756
2529436
240
501463
1030196
295292
2038472
360
449943
856206
408216
1636536In a biological sense Model2 implies that the inhibited control of the Raf-MEK pathway is directed via ERK, rather than ERK and RKIP, proteins by phosphorylating MEK at T292 binding site or by phosphorylating SOS, as long as the EGR receptor triggers the activation of the pathway. Moreover, the chains of reactions which produce, particularly, Shc-Grb2-SOSm, Grb2-SOS, ERK.p2, and Ras.GTP are very quick in the sense that Raf.A-Ras.GTPm, ERK.p1, and Shcm proteins, which are used in the production of those proteins, cannot stay longer in the cell to count them as separate species. On the other hand although the knowledge about the MAPK system claims the existence of all kinds of RKIP proteins in the system [32–34], Model 2 describes the system without any RKIP species.Once the selection of the best-fitted model we compare it with the proposal models of the system in [35] which are suggested by [36–38]. The main differences among these models are the function of the EGF receptor during the activation, the existence of Shc substrates and negative-feedback loops. Among these alternative models, Model 2 supports the model of [36]. The main feature of that model is the existence of the Ras activation with Shc and without Shc species, the degradation of EGRF, and core cascade of ERK via Raf and MEK. Moreover, different from this model, Model 2 includes the negative feedback phosphorylation of SOS by active ERK and defines the system with less species and reactions. Schoeberl et al. [36] model describes the pathway with 94 reactions and 125 reactions. On the other hand to the best of our knowledge since none of this model is validated by a real time-course dataset within a stochastic modelling, there have not yet any wet-data/approaches which can verify Model 2 or other suggested models.On the other hand, if we only evaluate the estimated values of reaction rates of Model2, we see that the dissociation type of reactions is the fast reactions and among other dissociations, the ones within the inactive Raf and active Ras complex (reaction 15 in Model 2) as well as within the active ERK, RSK, and transcription factor (reaction 39 in Model 2) have the highest speed. Whereas the reactions which regulate the Raf-MEK cascade by the inhibition of ERK (reactions 20 and 23 in Model 2) occur slowly. Figure 2 is shown as example from the posterior distributions of some reaction rates and the trace plots of their MCMC outputs. Figure 3, on the one hand, briefly describes Model 2. In this figure the thickness of the interaction indicates the speed of the associated reactions.Examples of the probability distribution (a, b, and c) and trace plots (d, e, and f) of reaction rates 2 (which indicates the simultaneous recruitment of Shc and Grb2-SOS complex from the cytosol to the cell membrane by the recruitment of EGFR (a and d)), 16 (which shows the activation of MEK proteins by active Raf (b and e)), and 38 (which refers to the dissociation of active ERK and RSK complex (c and f)) of Model2 under Δt=5 after 85,000 burn-in runs.
(a)
(b)
(c)
(d)
(e)
(f)Figure 3
Simple representation of the MAPK/ERK pathway via Model2 under Δt=5. The thickness of the arrow indicates the speed of the interaction. The thick arrows represent fast interactions (estimated reaction rates have power 10-1), only simple arrows describe moderately slow reactions (estimated reaction rates have power 10-2 or 10-3), and finally dotted lines show very slow reactions among substrates (estimated reaction rates have power 10-4 or 10-5).Figure 4
Simple representation of the structure of the MAPK/ERK pathway taken from Figure 2 in the study of [21].Figure 5
Plot of the western blotting data of the MAPK pathway observed at 8 time points according to the given molecular populations in Table5.On the other hand, [39] shows that the single element update of states causes poor mixing due to the high correlation within the latent data. Thereby as an extension of this study, we can investigate the possibility of block updates of missing states [40, 41] to get better mixing. We consider that this suggested approach improves the current approximation technique in the class of stochastic methods where the optimal solution of the system cannot be found globally. Moreover, as shown in Figure 2, the curves may also imply fractal time series [42, 43] in which the system can be also modelled via the fractional Brownian motion (FBM). In FBM, the density of the time series has heavy-tailed distribution and indicates slowly decreasing autocorrelation functions with power spectrum of power-law type that can be observed as different cell signalling processes [44] as well as earthquake modelling and financial time series analysis. The power-law spectrum is seen if the process shows the self-similarity, that is, a time segment in the process presents the same behaviour as any segments of other time scales, and stationarity of its increments in diffusion model [45, 46]. If the self-similarity is found locally, rather than globally as in FBM, the system can be modelled by a more general modelling approach, called the multifractional Brownian motion [45]. On the other side if we deal with the change in increment of the Gaussian process, the system can be described by the fractional Gaussian model (FGM) [47–50]. Hereby as the extension of the diffusion modelling, the high correlation in the system can be solved by both FBM and FGM models that are specifically designed for such short- or long-correlated time series datasets.


## 5. Conclusion

We have described the MAPK pathway by quasireactions and applied multiple parameterizations for the same protein in different localizations and different phosphorylations. This reaction set has been used for estimating reaction rates of the real time-course data under the assumption that the observed substrates possess the measurement errors. In inference of the model parameters we have used the nonlinear Euler approximation combined with a data augmentation step and implemented it within a Metropolis-within-Gibbs algorithm by taking into account the singularity of the system at different stages of the estimation. This method is one of the best-known mathematical methods in system and computational biology. The alternative models for the western blot MAPK data were compared via BIC and DIC techniques. Both results have shown Model2 as the best-selected model which fits the available experimental data. In this study even though we have been interested in the MAPK pathway as a result of its leading role in the cellular life cycle, the methods which we have developed for this system can be applied in the analysis of any large biochemical systems with many genes used in the biochemical and bioengineering researches.



---

*Source: 752631-2012-11-28.xml*