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.

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 [36] 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 notation . 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); (b)(c)(d);

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.) 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 ’s () which are the parameters to be estimated.

For inferring , 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), (b), (c), (d) Models 1 and 2,(e) Models 3, 4, and 5,

where Ras.GDP and Ras.GTP represent inactive and active Ras substrates, respectively. and denote MEK proteins phosphorylated by PAK and ERK proteins, in the given order. Raf.I and designate the inactive Raf in the cytosol and near the cell membrane, respectively, and finally Raf 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, MEK, , and 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 , where and 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 BG.

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 equations are summarized as resulting in the cancelation of 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, MEK.p2, ERK, and ERK.p2) and 14 unobserved substrates.

Model 2. It is the extension of Model 1 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. Model 1 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, 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 Model 3 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 molecules in each species at time , , is continuous, the stochastic model can be converted to a differential equation model. Indeed, can be further expanded via a Taylor series expansion and the change in states of each species at 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: in which and are mean, or drift, and variance, or diffusion, matrices, respectively, both depending on the state at time and the parameter vector explicitly. is the total number of substrates and states the total number of reactions in the system. is an -dimensional independent identically distributed Brownian random vector over time . is the net effect matrix which is the difference between the stoichiometric coefficients of products and reactants of each reaction. Finally, 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 matrix is composed of both observed species and unobserved or missing species at given time point , hereby where each at consists of totally 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 , 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.

Finally, in our model we assume that each observation at time 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 under measurement error, , via , where denotes the measurement error coming from normal distribution with mean and variance, also called the tuning parameter, for all time and observed substrates, that is, . In the sampling, since the tuning parameter has an effect on the acceptance rate of , 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 can be found in [21].

In inference of the parameters , the positivity of rate constants () is maintained by independent exponential priors with rate 1. For the prior of the measurement error , 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 . In our block scheme we divide the -dimensional vector of into small and equally sized groups with a dimension and then simulate each -dimensional group sequentially. The size of the block update dimension in each model is varied. We set 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 .

When we renew (), we begin to update column by column whose entries are composed of partially observed and completely augmented data. The candidate values of each 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 from a normal distribution with mean and variance via · denotes the observed values of the substrates at the given time 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 , 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, namely structural and incidental dependence. The former is originated from the linear dependence of the substrates in the net effect matrix . 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 or any numeric problems resulting from where is the -dimensional hazard vector of the system for a given and , and represents the diagonal matrix of . To work with the incidental dependency we check whether a candidate value for either missing states of 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 the western 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 Appendices BF. 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 and . Using PAM clustering [30] with 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 , 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 ’s, it is observed that the increase of the augmented states by going from to 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 than in . Moreover we get better mixing in some of the estimates under 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 min (12 latent states from to min and 24 latent states from to min), hereby, the dependency between the parameters and those states increases. On the other side, since the number of augmented states is not large from to min, 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 to min). 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.

As 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 (Section 2) 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 Table 1 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 than the estimates under , we still choose the results from for the model selection seeing that the estimates of the associated measurement errors are lower that the ones under . 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.

In a biological sense Model 2 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-, Grb2-SOS, ERK.p2, and Ras.GTP are very quick in the sense that Raf.A-Ras., ERK.p1, and 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 [3234], 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 [3638]. 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 Model 2, 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.

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) [4750]. 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 Model 2 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.

Appendices

A. Derivation of Candidate Generators for Partially Observed States

In order to estimate the stochastic reaction rate constants of the MAPK/ERK pathway (Figure 4) when the states of the system are partially observed, we use the candidate generators whose details of the derivations are presented as follows.

Under the assumption that the observed values involves the measurement error and generated by where , the mean of conditional on , denoted as , is derived via from the multivariate normal theories. Here , , and . Moreover, and . In that equation, refers to the correlation between and and equals to in which and . Furthermore, and , while () describes the transpose of the given vector. By substituting into (A.1), we get where has full rank [51].

Since in our estimation is similar to the idea of , , and , hereby (A.2) can be written as

On the other hand, in order to derive the conditional variance of , , under , we use from the consequence of the multivariate normal theory. In (A.4) seeing that , , and , the final is found via in which and . , similar to , has full rank [51].

Therefore, from , the candidate generator of conditional on , , which converges pointwise to when is sampled from with mean and variance . The candidate which is combination of both proposal and , sequentially, that is, , is accepted according to the acceptance probability

For the estimation of the stochastic reaction rate constants, we use the western blotting data presented in Table 5 and plotted in Figure 5. The data are gathered from the School for Cancer Studies in the University of Glasgow and are already published in the study of [20]. In inference of the parameters in all suggested MAPK/ERK models, we arbitrarily equate the lowest intensity to 10 molecules and extrapolate the remaining proportionally.

We use the following list of reactions for estimating the reaction rates of Model 1. The estimated rates and associated statistics are summarized in Table 7. The list of substrates, on the other hand, is given in Table 6.(1)Grb2 + SOS Grb2-SOS (2)EGFR + Shc + Grb2-SOS EGFR + Shc-Grb2-(3)EGFR + Grb2-SOS EGFR + Grb2-(4)Shc-Grb2- Shc + Grb2-SOS (5)Grb2-SOS Grb2 + SOS (6)Shc-Grb2-  + Ras.GDP Shc-Grb2- Ras.GTP (7)Grb2- Ras.GDP Grb2- Ras.GTP (8)Ras.GTP Ras.GDP (9)GAP + Ras.GTP GAP + Ras.GDP (10)Raf + PP2A Raf.I + PP2A (11)Raf.I + Ras.GTP Raf. Ras.GTP (12)Raf. Ras.GTP + PAK Raf. Ras.GTP + PAK (13)PP5 + Raf. PP5 + Raf.(14)Raf. Raf (15)Raf.I-Ras. Raf. Ras.GTP (16)Raf. MEK Raf. MEK.p2 (17)PAK + MEK PAK (18) Raf. MEK.p2 + Raf.(19)MEK.p2 + ERK MEK.p2 + ERK.p2 (20)ERK.p2 + MEK ERK.p2 (21) Raf. MEK.p2 + Raf.(22)ERK.p2 + Shc-Grb2- ERK.p2 + Shc- SOS (23)ERK.p2 + Grb2- ERK.p2 +    + SOS (24)Shc- Shc + Grb2 (25) Grb2 (26)ERK.p2 + TF ERK.p2-TF.p2 (27)ERK.p2-TF.p2 + c-Fos.DNA ERK.p2-TF.p2 + c-Fos.DNA + c-Fos.RNA + c-Fos (28)ERK.p2 + c-Fos ERK.p2 + c-Fos.p (29)ERK.p2-TF.p2 + MKP.DNA ERK.p2-TF.p2 + MKP.DNA + MKP.RNA + MKP (30)MKP + ERK.p2 MKP + ERK (31)ERK.p2 + RSK ERK.p2-RSK.A (32)ERK.p2-RSK.A + TF ERK.p2-RSK.A-TF.p2 (33)ERK.p2-RSK.A-TF.p2 + c-Fos.DNA ERK.p2-RSK.A-TF.p2 + c-Fos + c-Fos.DNA + c-Fos.RNA (34)ERK.p2-RSK.A + c-Fos ERK.p2-RSK.A + c-Fos.p (35)ERK.p2-RSK.A-TF.p2 + MKP.DNA ERK.p2-RSK.A-TF.p2 + MKP + MKP.DNA + MKP.RNA (36)MKP + ERK.p2-RSK.A MKP + ERK + RSK (37)ERK.p2-TF.p2 ERK.p2 + TF (38)ERK.p2-RSK.A ERK.p2 + RSK (39)ERK.p2-RSK.A-TF.p2 ERK.p2-RSK.A + TF

In Model 2, basically, we use the same reaction list given for Model 1. But different from that model, we also add the reaction of the EGFR degradation which is denoted by as the reaction. The estimated rate constants concerning Model 2 and other statistics are summarized in Table 8. Apart from the EGF receptor, the list of substrates and their classes are the same as given in Table 6. Under this model the EGF receptor is counted as independent, rather than dependent, term.

We apply the following list of reactions to infer the reaction rates of Model 3. The estimated reaction rates and calculated statistics are presented in Table 9. Table 10 gives the list of substrates used in inference.(1)Grb2 + SOS Grb2-SOS (2)EGFR + Shc EGFR +  (3)EGFR + Grb2-SOS EGFR + Grb2-(4)   + Grb2- Shc-Grb2-(5)Shc-Grb2-  + Grb2-(6) Shc (7)Grb2- Grb2-SOS (8)Grb2-SOS Grb2 + SOS (9)Shc-Grb2-  + Ras.GDP Shc-Grb2- Ras.GTP (10)Grb2- Ras.GDP Grb2- Ras.GTP (11)Ras.GTP Ras.GDP (12)GAP + Ras.GTP GAP + Ras.GDP (13)Raf + PP2A Raf.I + PP2A (14)Raf.I + Ras.GTP Raf.  + Ras.GTP (15)Raf.  + Ras.GTP Raf.I-Ras.(16)Raf.I-Ras. PAK Raf.A-Ras. PAK (17)Raf.A-Ras. Raf.  + Ras.GTP (18)PP5 + Raf. PP5 + Raf.(19)Raf. Raf (20)Raf.I-Ras. Raf. Ras.GTP (21)Raf. MEK Raf.   + MEK.p2 (22)PAK + MEK PAK (23)  + Raf. MEK.p2 + Raf.(24)MEK.p2 + ERK MEK.p2 + ERK.p1 (25)MEK.p2 + ERK.p1 MEK.p2 + ERK.p2 (26)ERK.p2 + MEK ERK.p2 +  (27)  + Raf. MEK.p2 + Raf.(28)ERK.p2 + Shc-Grb2- ERK.p2 + Shc- SOS (29)ERK.p2 + Grb2- ERK.p2 +    + SOS (30)Shc- Shc + Grb2 (31) Grb2 (32)ERK.p2 + TF ERK.p2-TF.p2 (33)ERK.p2-TF.p2 + c-Fos.DNA ERK.p2-TF.p2 + c-Fos.DNA + c-Fos.RNA (34)c-Fos.RNA c-Fos (35)ERK.p2 + c-Fos ERK.p2 + c-Fos.p (36)ERK.p2-TF.p2 + MKP.DNA ERK.p2-TF.p2 + MKP.DNA + MKP.RNA (37)MKP.DNA MKP (38)MKP + ERK.p2 MKP + ERK (39)ERK.p2 + RSK ERK.p2-RSK.A (40)ERK.p2-RSK.A + TF ERK.p2-RSK.A-TF.p2 (41)ERK.p2-RSK.A-TF.p2 + c-Fos.DNA ERK.p2-RSK.A-TF.p2 + c-Fos.DNA + c-Fos.RNA (42)ERK.p2-RSK.A + c-Fos ERK.p2-RSK.A + c-Fos.p (43)ERK.p2-RSK.A-TF.p2 + MKP.DNA ERK.p2-RSK.A-TF.p2 + MKP.DNA + MKP.RNA (44)MKP + ERK.p2-RSK.A MKP + ERK + RSK (45)ERK.p2-TF.p2 ERK.p2 + TF (46)ERK.p2-RSK.A ERK.p2 + RSK (47)ERK.p2-RSK.A-TF.p2 ERK.p2-RSK.A + TF

In Model 4, apart from the additional reaction of the EGFR degradation as the reaction, we apply the same list of reactions presented for Model 3. Similarly we use the same list of substrates stated for Model 3. But in this case the receptor of the EGF protein (EGFR) is grouped with the independent, rather than dependent, species as found in Model 2. The estimates of rate constants and corresponding statistics are given in Table 11.

We take the following list of reactions to infer the rate constants of Model 5. Apart from those reactions of activations, we list all associated degradations of the proteins after their dissociations in Table 12. Table 14 summarizes the estimated model parameters merely via the following reactions list. On the other hand, the substrates used in inference are listed in Table 13.(1)Grb2 + SOS Grb2-SOS (2)EGFR + Shc EGFR +  (3)EGFR + Grb2-SOS EGFR + Grb2-(4)  + Grb2- Shc-Grb2-(5)Shc-Grb2-  + Grb2-(6) Shc (7)Grb2- Grb2-SOS (8)Grb2-SOS Grb2 + SOS (9)Shc-Grb2-  + Ras.GDP Shc-Grb2-  + Ras.GTP (10)Grb2- Ras.GDP Grb2-  + Ras.GTP (11)Ras.GTP Ras.GDP (12)GAP + Ras.GTP GAP + Ras.GDP (13)Raf + PP2A Raf.I + PP2A (14)Raf.I + Ras.GTP Raf.  + Ras.GTP (15)Raf.  + Ras.GTP Raf.I-Ras.(16)Raf.I-Ras.  + PAK Raf.A-Ras.  + PAK (17)Raf.A-Ras. Raf. Ras.GTP (18)PP5 + Raf. PP5 + Raf.(19)Raf. Raf (20)Raf.I-Ras. Raf.  + Ras.GTP (21)Raf.  + MEK Raf.  + MEK.p2 (22)PAK + MEK PAK +  (23)  + Raf. MEK.p2 + Raf.(24)Raf.I + RKIP Raf.I-RKIP (25)Raf.I-RKIP + Ras.GTP Raf.I-  + Ras.GTP (26)Raf.I-  + Ras.GTP Raf.I-RKIP-Ras.(27)MEK + RKIP MEK-RKIP (28)  + RKIP -RKIP (29)+ RKIP -RKIP (30)MEK.p2 + RKIP MEK.p2-RKIP (31)PKC + Raf.I-RKIP PKC + Raf.I + RKIP.p (32)ERK.p2 + Raf.I-RKIP ERK.p2 + Raf.I + RKIP.p (33)Raf.I-RKIP-Ras. Raf.I-  + Ras.GTP (34)Raf.I- Raf.I-RKIP (35)MEK-RKIP MEK + RKIP (36)-RKIP   + RKIP (37)-RKIP   + RKIP (38)MEK.p2-RKIP MEK.p2 + RKIP (39)RKIP.p RKIP (40)MEK.p2 + ERK MEK.p2 + ERK.p1 (41)MEK.p2 + ERK.p1 MEK.p2 + ERK.p2 (42)MEK.p2-RKIP + ERK MEK.p2-RKIP + ERK.p1 (43)MEK.p2-RKIP + ERK.p1 MEK.p2-RKIP + ERK.p2 (44)ERK.p2 + MEK ERK.p2 +  (45)  + Raf. MEK.p2 + Raf.(46)ERK.p2 + Shc-Grb2- ERK.p2 + Shc-  + SOS (47)ERK.p2 + Grb2- ERK.p2 +    + SOS (48)Shc- Shc + Grb2 (49) Grb2 (50)ERK.p2 + TF ERK.p2-TF.p2 (51)ERK.p2-TF.p2 + c-Fos.DNA ERK.p2-TF.p2 + c-Fos.DNA + c-Fos.RNA (52)c-Fos.RNA c-Fos (53)ERK.p2 + c-Fos ERK.p2 + c-Fos.p (54)ERK.p2-TF.p2 + MKP.DNA ERK.p2-TF.p2 + MKP.DNA + MKP.RNA (55)MKP.DNA MKP (56)MKP + ERK.p2 MKP + ERK (57)ERK.p2 + RSK ERK.p2-RSK.A (58)ERK.p2-RSK.A + TF ERK.p2-RSK.A-TF.p2 (59)ERK.p2-RSK.A-TF.p2 + c-Fos.DNA ERK.p2-RSK.A-TF.p2 + c-Fos.DNA + c-Fos.RNA (60)ERK.p2-RSK.A + c-Fos ERK.p2-RSK.A + c-Fos.p (61)ERK.p2-RSK.A-TF.p2 + MKP.DNA ERK.p2-RSK.A-TF.p2 + MKP.DNA + MKP.RNA (62)MKP + ERK.p2-RSK.A MKP + ERK + RSK (63)ERK.p2-TF.p2 ERK.p2 + TF (64)ERK.p2-RSK.A ERK.p2 + RSK (65)ERK.p2-RSK.A-TF.p2 ERK.p2-RSK.A + TF (66)EGFR

G. List of Protein States Used in the MAPK Pathway

We use the following substrates in the description of the MAPK pathway with different models.(1) Ras protein states.(a)Ras.GDP: the inactive Ras protein near the cell membrane.(b)Ras.GTP: the active Ras near the cell membrane.(2)Raf protein states(a)Raf: the inactive and nonphosphorylated Raf protein in the cytosol.(b)Raf.: the active Raf phosphorylated on the S338 and the S471 binding sites near the cell membrane.(c)Raf.A-Ras.: the complex of the active Raf and Ras.GTP near the cell membrane.(d)Raf.I: the inactive Raf phosphorylated on the S259 binding site in the cytosol.(e)Raf.: the inactive Raf phosphorylated on the S259 binding site and recruited from the cytosol to the cell membrane by the Ras.GTP protein.(f)Raf.I-Ras.: the complex of the inactive Raf and Ras.GTP near the cell membrane.(g)Raf.I-RKIP: the complex of the inactive Raf and RKIP, whose binding site is S338, in the cytosol.(h)Raf.I-: the complex of the inactive Raf and RKIP which is recruited to the membrane by the Ras.GTP protein.(i)Raf.I-RKIP-Ras.: the complex of the inactive Raf, RKIP, and Ras.GTP near the cell membrane.(3)MEK protein states.(a)MEK: the inactive and nonphosphorylated MEK protein in the cytosol.(b): the inactive MEK in the cytosol which is monophosphorylated by the activator PAK on the S298 binding site.(c): the inactive MEK in the cytosol which is monophosphorylated by the active ERK on the T292 binding site.(d)MEK.p2: the double-phosphorylated MEK (active MEK) on the S218 and S222 binding sites in the cytosol.(e)MEK-RKIP: the complex of the MEK and RKIP proteins in the cytosol.(f)-RKIP: the complex of the and RKIP proteins in the cytosol.(g)-RKIP: the complex of the and RKIP proteins in the cytosol.(h)MEK.p2-RKIP: the complex of the active MEK.p2 and RKIP proteins in the cytosol.(4)ERK protein states.(a)ERK: the inactive and nonphosphorylated ERK protein in the cytosol.(b)ERK.p1: the inactive, monophosphorylated ERK in the cytosol.(c)ERK.p2: the double-phosphorylated ERK (active ERK) in the cytosol.(d)ERK.p2-RSK.A: the complex of the active ERK and active RSK, which is activated by the ERK protein, in the nucleus.(e)ERK.p2-RSK.A-TF.p2: the complex of the active ERK, active RSK, and double-phosphorylated transcription factor in the nucleus.(f)ERK.p2-TF.p2: the complex of the active ERK and a transcription factor (like Elk or SAP proteins), which is double-phosphorylated by the active ERK, in the nucleus.(5)Grb2, Shc, and SOS protein states.(a)Grb2: A protein in the cytosol.(b): the Grb2 protein near the cell membrane after the dissociation of the SOS protein by the active ERK.(c)Grb2-SOS: the complex of the Grb2 and SOS proteins in the cytosol.(d)Grb2-: the complex of the Grb2 and SOS proteins near the cell membrane, where it is able to activate the Ras protein.(e)Shc: a protein in the cytosol.(f): the Shc protein near the cell membrane after the activation of the EGF receptor.(g)Shc-: the complex of the Shc and Grb2 proteins near the cell membrane after the dissociation of the SOS protein by the active ERK.(h)Shc-Grb2-: the complex of the Shc, Grb2, and SOS proteins near the cell membrane, where it is able to activate the Ras protein.(i)SOS: a protein, which is an exchange factor, in the cytosol.(6)c-Fos and MKP protein states.(a)c-Fos: a protein in the nucleus.(b)c-Fos.DNA: the gene sequence of the c-Fos protein.(c)c-Fos.p: the c-Fos gene phosphorylated by the ERK protein.(d)c-Fos.RNA: the transcription of the c-Fos gene into the messenger RNA (mRNA).(e)MKP: a protein in the cytosol.(f)MKP.DNA: the gene sequence of the MKP protein.(g)MKP.RNA: the transcription of the MKP gene into mRNA.(7)Other proteins(a)EGF: a protein which triggers the activation of the pathway by attaching its receptor (EGFR) in the cell membrane.(b)EGFR: a receptor that is equated with activated tyrosine kinase receptors. (c)GAP: a protein near the cell membrane.(d)PAK: a protein near the cell membrane.(e)PKC: a protein in cytosol.(f)PP2A: a protein near the cell membrane or in cytosol.(g)PP5: a protein near the cell membrane.(h)RKIP: a protein in the cytosol.(i)RKIP.p: the RKIP monophosphorylated either by the PKC or ERK protein on the binding sites S153 and S99, respectively.(j)RSK: an inactive protein in the cytosol.(k)RSK.A: the active RSK, which is activated by the ERK.p2 protein.(l)TF: a transcription factor (like Elk or SAP proteins), which will be double-phosphorylated by the active ERK, in the nucleus.

Acknowledgments

The authors would like to thank Professor Walter Kolch from the Conway Institute of Biomolecular and Biomedical Research for his valuable explanation about the MAPK/ERK pathway and providing the data. Also they would like to thank the GRID-TR service as they provided them with use the high performance computer.