www.prismmodelchecker.org

MAPK Cascade

(Huang and Ferrell)

Contents

Related publications: [KNP08a]

Introduction

The MAP (Mitogen-Activated Protien) Kinases are involved in a pathway through which information is sent to the nucleus. It is one of the most important signalling pathways since it plays a pivotal role in the molecular-signalling that governs the growth, proliferation and survival of many cell types. The MAPK cascade consists of a MAPK Kinase Kinase (MAPKKK), MAPK Kinase (MAPKK) and a MAPK. The cascade is initialised through the phosphorylation of MAPKKK, which then activates MAPKK through the phosphorylation at two serine residues and this then activates MAPK through the phosphorylation at theronine and tyrosine residues. The initialisation can be caused by a diverse set of stimuli including growth factors, neurotransmitters and cytokines.

The Pathway

The model of the pathway we consider is taken from [HF96]. Its structure is given below.

image of pathway

Following [HF96], we assume that the phosphorylation of both MAPK and MAPKK occur in two distributed steps. For example, when MAPK collides with its activator (MAPKK-PP) the first phosphorylation (MAPK-P) occurs and the activator is released. The phosphorylated MAPK must then collide again with its activator for the second phosphorylation (MAPK-PP) to occur. The deactivation of phosphorylated MAPK and MAPKK is caused by the corresponding phosphatase, while the activation and deactivation of MAPKKK is through the enzymes E1 and E2 respectively. To simplify the presentation we denote MAPK, MAPKK and MAPKKK by K, KK and KKK respectively.

1. MAPKKK is activated through enzyme E1
         KKK + E1 → KKK:E1 (a1 = 1 nM-1s-1)
         KKK + E1 ← KKK:E1 (d1 = 150 s-1)
         KKK:E1 → KKK* + E1 (k1 = 150 s-1)
2. MAPKKK is deactivated through enzyme E2
         KKK* + E2 → KKK:E2 (a2 = 1 nM-1s-1)
         KKK* + E2 ← KKK:E2 (d2 = 150 s-1)
         KKK*:E2 → KKK + E2 (k2 = 150 s-1)
3. MAPKK is activated by MAPKKK*
         KK + KKK* → KK:KKK* (a3 = 1 nM-1s-1)
         KK + KKK* ← KK:KKK* (d3 = 150 s-1)
         KK:KKK* → KK-P + KKK* (k3 = 150 s-1)
4. MAPKK-P is deactivated by MAPKK phosphatase
         KK-P + KK-Ptase → KK-P:KK-Ptase (a4 = 1 nM-1s-1)
         KK-P + KK-Ptase ← KK-P:KK-Ptase (d4 = 150 s-1)
         KK-P:KK-Ptase → KK + KK-Ptase (k4 = 150 s-1)
5. MAPKK-P is activated by MAPKKK*
         KK-P + KKK* → KK-P:KKK* (a5 = 1 nM-1s-1)
         KK-P + KKK* ← KK-P:KKK* (d5 = 150 s-1)
         KK-P:KKK* → KK-PP + KKK* (k5 = 150 s-1)
6. MAPKK-PP is deactivated by MAPKK phosphatase
         KK-PP + KK-Ptase → KK-PP:KK-Ptase (a6 = 1 nM-1s-1)
         KK-PP + KK-Ptase ← KK-PP:KK-Ptase (d6 = 150 s-1)
         KK-PP:KK-Ptase → KK-P + KK-Ptase (k6 = 150 s-1)
7. MAPK is activated by MAPKK-PP
         K + KK-PP → K:KK-PP (a7 = 1 nM-1s-1)
         K + KK-PP ← K:KK-PP (d7 = 150 s-1)
         K:KK-PP → K-P + KK-PP (k7 = 150 s-1)
8. MAPK-P is deactivated by MAPK phosphatase
         K-P + K-Ptase → K-P:K-Ptase (a8 = 1 nM-1s-1)
         K-P + K-Ptase ← K-P:K-Ptase (d8 = 150 s-1)
         K-P:K-Ptase → K + K-Ptase (k8 = 150 s-1)
9. MAPK-P is activated by MAPKK-PP
         K-P + KK-PP → K-P:KK-PP (a9 = 1 nM-1s-1)
         K-P + KK-PP ← K-P:KK-PP (d9 = 150 s-1)
         K-P:KK-PP → K-PP + KK-PP (k9 = 150 s-1)
10. MAPK-PP is deactivated by MAPK phosphatase
         K-PP + K-Ptase → K-PP:K-Ptase (a10 = 1 nM-1s-1)
         K-PP + K-Ptase ← K-PP:K-Ptase (d10 = 150 s-1)
         K-PP:K-Ptase → K-P + K-Ptase (k10 = 150 s-1)

The kinetic rates are based on the data presented in [HF96] where it is assumed the Km values (Km = (dm + km) / am) for the phosphorylation and dephosphorylation of MAPK, MAPKK and MAPKK all equal 300 nM.

The PRISM Model

When building a PRISM model of a biological pathway, it is possible to construct an individual-based model which provides a detailed model of the evolution of individual molecular components. However, taking this approach comes at a cost: it will inevitably suffer from the well known state-space explosion problem where, as the complexity of the system increases, the state space of the underlying model grows exponentially.

An alternative is to employ a population-based approach where the number of each type of molecule or species is modelled, rather than the state of each individual component. Such an approach leads to a much smaller state-space (see for example [HKN+08]) while still including sufficient detail to express the properties of interest. For these reasons, it is this approach that we use here.

We now describe the specification in PRISM of the MAPK cascade from the previous section. Each of the basic elements of the pathway (MAPKKK, MAPKK, MAPK, MAPKK Phosphatase, MAPK Phosphatase, E1 and E2), is represented by a separate PRISM module and synchronisation between modules is used to model reactions which involve interactions between multiple elements. However, the different forms which each can take (for example, which other compounds it is bound to) are represented by one or more variables within the module.

The PRISM model for the pathway is given below.

// MAPK cascade 
// taken from: C.-Y. Huang and J. Ferrell
// Ultrasensitivity in the mitgen-activated protein kinase cascade
// Proc. Natl. Acad. Sci. 93:10078-13100, 1996
// gxn/dxp 19/12/07

ctmc

// constants (number of elements of each species)
const int E1=1; // initial amount of E1
const int E2=1; // initial amount of E2
const int M=1;  // initial amount of MAPK phosphatase and MAPKK phosphatase
const int N;  // initial amount of MAPK/MAPKK and MAPKKK

// reaction rates (suppose volume proportional to the number of elements of MAPK)
const double a1=1/N;
const double d1=150;
const double k1=150;
const double a2=1/N;
const double d2=150;
const double k2=150;
const double a3=1/N;
const double d3=150;
const double k3=150;
const double a4=1/N;
const double d4=150;
const double k4=150;
const double a5=1/N;
const double d5=150;
const double k5=150;
const double a6=1/N;
const double d6=150;
const double k6=150;
const double a7=1/N;
const double d7=150;
const double k7=150;
const double a8=1/N;
const double d8=150;
const double k8=150;
const double a9=1/N;
const double d9=150;
const double k9=150;
const double a10=1/N;
const double d10=150;
const double k10=150;

//------------------------------------------------------------------------------
// enzymes

module E1
	
	e1 : [0..E1] init E1; // amount of enzyme E1
	
	// reaction 1 (MAPKKK is activated by E1)
	[a_kkk_e1] e1>0  -> e1 : (e1'=e1-1);
	[d_kkk_e1] e1<E1 ->  1 : (e1'=e1+1);
	[k_kkk_e1] e1<E1 ->  1 : (e1'=e1+1);
	
endmodule

// construct E2 by renaming E1
module E2 = E1[
	e1=e2,
	E1=E2,
	a_kkk_e1=a_kkk_e2,
	d_kkk_e1=d_kkk_e2,
	k_kkk_e1=k_kkk_e2 ]
endmodule

//------------------------------------------------------------------------------
// phosphatases

// mapk phosphatases
module KPTASE
	
	kptase : [0..M] init M; // amount of MAPK Phosphatase
	
	// reactions 8 and 10 (MAPK/MAPK-P is deactivated by MAPK Phosphatase)
	[a_k_ptase] kptase>0 -> kptase : (kptase'=kptase-1);
	[d_k_ptase] kptase<M -> 1 : (kptase'=kptase+1);
	[k_k_ptase] kptase<M -> 1 : (kptase'=kptase+1);
	
endmodule

// construct mapkk phosphatases by renaming module for mapk phosphatases
module KKPTASE = KPTASE[
	kptase= kkptase,
	a_k_ptase=a_kk_ptase,
	d_k_ptase=d_kk_ptase,
	k_k_ptase=k_kk_ptase ]
endmodule

//------------------------------------------------------------------------------
// mapks

// mapk
module MAPK
	
	k         : [0..N] init N; // quantity of MAPK
	k_kkpp    : [0..N] init 0; // quantity of MAPK:MAPKK-PP
	kp        : [0..N] init 0; // quantity of MAPK-P
	kp_kkpp   : [0..N] init 0; // quantity of MAPK-P:MAPKK-PP
	kp_ptase  : [0..N] init 0; // quantity of MAPK-P:MAPK Phosphatase
	kpp       : [0..N] init 0; // quantity of MAPK-PP
	kpp_ptase : [0..N] init 0; // quantity of MAPK-PP:MAPK Phosphatase
	
	// reaction 7 (MAPK is activated by MAPKK-PP)
	[a_k_kk] k>0 & k_kkpp<N  -> a7*k      : (k_kkpp'=k_kkpp+1) & (k'=k-1);
	[d_k_kk] k<N & k_kkpp>0  -> d7*k_kkpp : (k_kkpp'=k_kkpp-1) & (k'=k+1);
	[k_k_kk] k_kkpp>0 & kp<N -> k7*k_kkpp : (k_kkpp'=k_kkpp-1) & (kp'=kp+1);
	// reaction 8 (MAPK-P is deactivated by MAPK Phosphatase)
	[a_k_ptase] kp>0 & kp_ptase<N -> a8*kp       : (kp_ptase'=kp_ptase+1) & (kp'=kp-1);
	[d_k_ptase] kp<N & kp_ptase>0 -> d8*kp_ptase : (kp_ptase'=kp_ptase-1) & (kp'=kp+1);
	[k_k_ptase] kp_ptase>0 & k<N  -> k8*kp_ptase : (kp_ptase'=kp_ptase-1) & (k'=k+1);
	// reaction 9 (MAPK-P is activated by MAPKK-PP)
	[a_k_kk] kp>0 & kp_kkpp<N  -> a9*kp      : (kp_kkpp'=kp_kkpp+1) & (kp'=kp-1);
	[d_k_kk] kp<N & kp_kkpp>0  -> d9*kp_kkpp : (kp_kkpp'=kp_kkpp-1) & (kp'=kp+1);
	[k_k_kk] kp_kkpp>0 & kpp<N -> k9*kp_kkpp : (kp_kkpp'=kp_kkpp-1) & (kpp'=kpp+1);
	// reaction 10 (MAPK-P is deactivated by MAPK Phosphatase)
	[a_k_ptase] kpp>0 & kpp_ptase<N -> a10*kpp       : (kpp_ptase'=kpp_ptase+1) & (kpp'=kpp-1);
	[d_k_ptase] kpp<N & kpp_ptase>0 -> d10*kpp_ptase : (kpp_ptase'=kpp_ptase-1) & (kpp'=kpp+1);
	[k_k_ptase] kpp_ptase>0 & kp<N  -> k10*kpp_ptase : (kpp_ptase'=kpp_ptase-1) & (kp'=kp+1);
	
endmodule

// mapkk
module MAPKK
	
	kk         : [0..N] init N; // quantity of MAPKK
	kk_kkkp    : [0..N] init 0; // quantity of MAPKK:MAPKKK-P
	kkp        : [0..N] init 0; // quantity of MAPKK-P
	kkp_kkkp   : [0..N] init 0; // quantity of MAPKK-P:MAPKKK-P
	kkp_ptase  : [0..N] init 0; // quantity of MAPKK-P:MAPKK Phosphatase
	kkpp       : [0..N] init 0; // quantity of MAPKK-PP
	kkpp_ptase : [0..N] init 0; // quantity of MAPKK-PP:MAPKK Phosphatase
	
	// reaction 3 (MAPKK is activated by MAPKKK*)
	[a_kk_kkk] kk>0 & kk_kkkp<N  -> a3*kk      : (kk_kkkp'=kk_kkkp+1) & (kk'=kk-1);
	[d_kk_kkk] kk<N & kk_kkkp>0  -> d3*kk_kkkp : (kk_kkkp'=kk_kkkp-1) & (kk'=kk+1);
	[k_kk_kkk] kk_kkkp>0 & kkp<N -> k3*kk_kkkp : (kk_kkkp'=kk_kkkp-1) & (kkp'=kkp+1);
	// reaction 4 (MAPKK-P is deactivated by MAPKK Phosphatase)
	[a_kk_ptase] kkp>0 & kkp_ptase<N -> a4*kkp       : (kkp_ptase'=kkp_ptase+1) & (kkp'=kkp-1);
	[d_kk_ptase] kkp<N & kkp_ptase>0 -> d4*kkp_ptase : (kkp_ptase'=kkp_ptase-1) & (kkp'=kkp+1);
	[k_kk_ptase] kkp_ptase>0 & kk<N  -> k4*kkp_ptase : (kkp_ptase'=kkp_ptase-1) & (kk'=kk+1);
	// reaction 5 (MAPKK-P is activated by MAPKKK*)
	[a_kk_kkk] kkp>0 & kkp_kkkp<N  -> a5*kkp      : (kkp_kkkp'=kkp_kkkp+1) & (kkp'=kkp-1);
	[d_kk_kkk] kkp<N & kkp_kkkp>0  -> d5*kkp_kkkp : (kkp_kkkp'=kkp_kkkp-1) & (kkp'=kkp+1);
	[k_kk_kkk] kkp_kkkp>0 & kkpp<N -> k5*kkp_kkkp : (kkp_kkkp'=kkp_kkkp-1) & (kkpp'=kkpp+1);
	// reaction 6 (MAPKK-P is deactivated by MAPKK Phosphatase)
	[a_kk_ptase] kkpp>0 & kkpp_ptase<N -> a6*kkpp       : (kkpp_ptase'=kkpp_ptase+1) & (kkpp'=kkpp-1);
	[d_kk_ptase] kkpp<N & kkpp_ptase>0 -> d6*kkpp_ptase : (kkpp_ptase'=kkpp_ptase-1) & (kkpp'=kkpp+1);
	[k_kk_ptase] kkpp_ptase>0 & kkp<N  -> k6*kkpp_ptase : (kkpp_ptase'=kkpp_ptase-1) & (kkp'=kkp+1);
	// reactions 7 and 9 (MAPK/MAPK-P is activated by MAPKPP)
	[a_k_kk] kkpp>0 -> kkpp : (kkpp'=kkpp-1);
	[d_k_kk] kkpp<N -> 1 : (kkpp'=kkpp+1);
	[k_k_kk] kkpp<N -> 1 : (kkpp'=kkpp+1);
	
endmodule

// mapkkk
module MAPKKK
	
	kkk     : [0..N] init N; // quantity of MAPKKK
	kkk_e1  : [0..N] init 0; // quantity of MAPKKK:E1
	kkkp    : [0..N] init 0; // quantity of MAPKKK*
	kkkp_e2 : [0..N] init 0; // quantity of MAPKKK*:E2
	
	// reaction 1 (MAPKKK is activated by E1)
	[a_kkk_e1] kkk>0 & kkk_e1<N  -> a1*kkk    : (kkk_e1'=kkk_e1+1) & (kkk'=kkk-1);
	[d_kkk_e1] kkk<N & kkk_e1>0  -> d1*kkk_e1 : (kkk_e1'=kkk_e1-1) & (kkk'=kkk+1);
	[k_kkk_e1] kkk_e1>0 & kkkp<N -> k1*kkk_e1 : (kkk_e1'=kkk_e1-1) & (kkkp'=kkkp+1);
	// reaction 2 (MAPKKK* is deactivated by E2)
	[a_kkk_e2] kkkp>0 & kkkp_e2<N -> a2*kkkp    : (kkkp_e2'=kkkp_e2+1) & (kkkp'=kkkp-1);
	[d_kkk_e2] kkkp<N & kkkp_e2>0 -> d2*kkkp_e2 : (kkkp_e2'=kkkp_e2-1) & (kkkp'=kkkp+1);
	[k_kkk_e2] kkkp_e2>0 & kkk<N  -> k2*kkkp_e2 : (kkkp_e2'=kkkp_e2-1) & (kkk'=kkk+1);
	// reactions 3 and 5 (MAPKK/MAPKK-P is activated by MAPKKP)
	[a_kk_kkk] kkkp>0 -> kkkp : (kkkp'=kkkp-1);
	[d_kk_kkk] kkkp<N -> 1 : (kkkp'=kkkp+1);
	[k_kk_kkk] kkkp<N -> 1 : (kkkp'=kkkp+1);
	
endmodule

//------------------------------------------------------------------------------
// reward structures

rewards "activated" // activated mapk
	true : kpp;
endrewards

rewards "activated_squared" // activated mapk squared (used to calculate standard deviation)
	true : kpp*kpp;
endrewards

rewards "percentage" //percentage activated mapk
	true : 100*(kpp/N);
endrewards

rewards "reactions" // reactions between mapk and mapkk
	[a_k_kk] true : 1;
	[d_k_kk] true : 1;
	[k_k_kk] true : 1;
endrewards

rewards "time" // time
	true : 1;
endrewards
View: printable version          Download: mapk_cascade.sm

We have included five reward structures in the description. The first reward structure (''activated'') assigns a state reward equal to the amount of MAPK that is activated while the second reward structure (''activated_squared'') assigns a state reward equal to the square of the amount MAPK that is activated. Since, the standard deviation of a random variable X equals the square root of its variance which equals:

E(X2) - E(X)2,

together these rewards can be used to compute the standard deviation (and variance) of random variable for the amount of activated MAPK at a time instant. The third (''activated'') assigns a state reward equal to the percentage of MAPK that is activated. The forth reward structure (''reactions'') assigns a reward of 1 to all transitions which correspond to a reaction between MAPK and MAPKK and the last (''time'') simply assigns a state reward of 1 to all states in the model.

The table below shows statistics for the MTBDD which represents each model we have built, in terms of the number of states, transitions, the total number of nodes and the number of leaves (terminal nodes).

N States Transitions Nodes Leaves
1 118 468 1,976 3
2 2,172 13,608 18,095 6
3 18,292 144,630 41,72810
4 99,535 910,872111,56314
5 408,366 4,138,848185,92820
6 1,373,026 15,015,264285,10625
7 3,979,348 46,167,582414,83533
810,276,461125,012,862740,27639

The table below gives the times taken to construct the models. There are two stages. Firstly, the system description (in our module based language) is parsed and converted to the MTBDD representing it. Secondly, reachability is performed to identify non-reachable states and the MTBDD is filtered accordingly. Reachability is performed using a BDD based fixpoint algorithm. The number of iterations for this process is also given.

N Time (s) Iter.s
10.09319
20.39225
31.29734
47.77544
519.2954
641.6464
789.5474
8217.384

Model Checking

We analyse this pathway through the following properties:

  • the amount of activated MAPK at time t and the standard deviation of this random variable;
  • the expected percentage of activated MAPK at time t;
  • the expected number of reactions between MAPK and MAKK by time t;
  • the expected time until all MAPK are activated at the same time instant.

These property can be computed by the CSL through the following formulae:

const double t; // time bound

// the expected amount of activated MAPK at time t
R{"activated"}=?[ I=t ]
// the expected amount squared of activated MAPK at time t (used to calculate standard deviation)
R{"activated_squared"}=?[ I=t ]
// the expected percentage of activated MAPK at time t
R{"percentage"}=?[ I=t ]
//the expected number of reactions between MAPK and MAKK by time t
R{"reactions"}=?[ C<=t ]
// the expected time until all MAPK are activated at the same time instant
R{"time"}=?[F kpp=N ]
View: printable version          Download: mapk_cascade.csl
  • The amount of activated MAPK at time t and the standard deviation of this random variable.

    In the graph below, we have plotted this measure as t varies for N=4 and N=8. The graphs also include the standard deviation of the random variable for the amount of activated MAPK at time t, drawn as a pair of dotted lines

    • N=4

      graph plotting the expected amount of activated MAPK at time t (N=4)
    • N=8

      graph plotting the expected amount of activated MAPK at time t (N=8)

    For the purposes of comparison, we also show results for the expected amount of activated MAPK computed using PRISM's discrete-event simulation engine. These are generated using very small numbers of simulation runs (10 and 100). Smoother approximations for the plots above can be obtained with higher numbers of runs.

    • N=4

      graph plotting simulation results for the amount of activated MAPK at time t (N=4)
    • N=8

      graph plotting simulation results for the amount of activated MAPK at time t (N=8)

    In the graph below we have ploted how the standard deviation of the random variable for the amount of activated MAPK at time t changes as both t and N vary

    graph plotting the standard deviation of the rv for the amount activated MAPK at time t
  • The expected percentage of activated MAPK at time t.

    In the graph below, we have plotted this measure as both t and N vary.

    graph plotting the expected percentage of activated MAPK at time t
  • The expected number of reactions between MAPK and MAKK by time t.

    In the graph below, we have plotted this measure as both t and N vary.

    graph plotting the expected number of reactions between MAPK and MAKK by time t
  • the expected time until all MAPK are activated at the same time instant.

    In the graph below we present the results obtained for this property as both N and L vary.

    graph plotting the expected time until all MAPK are activated at the same time instant

Case Studies