	PathwayLab Research Edition

	Features
	Examples
		BioCircuit
		A MAPK Cascade
		Integral Feedback Control
		Cell Cycle Control in Xenopus Frog Eggs

	Buy Online
	Trial Version
	For More Information
	Life & Medical Sciences Solutions

Stanford University

	Ask about this page
	Print this page
	Give us feedback
	Sign up for the Wolfram Insider

Example: A MAPK Cascade

This model is taken from B. N. Kholodenko et al. 2002 ¹. The model can be launched in PathwayLab Research Edition from the Start menu, but its pathway diagram is also shown as plain graphics in the figure below. The file MAPKCascade.nb^* used below was generated by Export to > Mathematica in the PathwayLab menu.

[Graphics:HTMLFiles/mapk_1.gif]

A mitogen-activated protein kinase (MAPK) cascade

Load the package.

In[1]:=

Load the model and store it in a variable.

In[2]:=

Out[2]=

Here are the state variables of the pathway.

In[3]:=

Out[3]=

These are the rate variables of the pathway.

In[4]:=

Out[4]=

This lists mass balance equations.

In[5]:=

Out[6]//TraditionalForm=

Once the model equations are available in Mathematica, it is easy to analyze them in various ways. Here we will compute the stoichiometric matrix.

This gives the right-hand side of the mass balance equations.

In[7]:=

Out[7]=

{v5[t] - v6[t] + v7[t] - v8[t], v6[t] - v7[t], v10[t] - v11[t], -v1[t] + v4[t], v1[t] - v2[t] + v3[t] - v4[t], v2[t] - v3[t], -v5[t] + v8[t], -v10[t] + v11[t] - v12[t] + v9[t], v12[t] - v9[t]}

This gives the stoichiometric matrix (in Mathematica 5.1 or later this can be directly accomplished by D[rhs,{#[t]&/@rates}]).

In[8]:=

Out[9]//MatrixForm=

Here are the mass balance equations in matrix notation.

In[10]:=

Out[10]=

Linear dependent rows of the stoichiometric matrix imply that there are conserved moieties in the system. In this particular case it is easy to find these linear combinations by inspection. The corresponding linear relations between the state derivatives show that certain combinations of states are constant over time.

Here are the linear combinations and a test that they map zero.

In[11]:=

Out[14]//MatrixForm=

These are the conservation equations.

In[15]:=

$qeqn = Thread[{MKK_tot, MAPK_tot, MKKK_tot} == {l1, l2, l3} . (#[t] &/@states)]$

Out[15]=

${MKK_tot == MKK[t] + MKK—P[t] + MKK—PP[t], MAPK_tot == MAPK[t] + MAPK—P[t] + MAPK—PP[t], MKKK_tot == MKKK[t] + MKKK—P[t] + MKKK—PP[t]}$

This gives a list of rules for the total amounts of MKKK, MKK, and MAPK.

In[16]:=

Out[16]=

${MKK_tot→180, MAPK_tot→360, MKKK_tot→200}$

By utilizing the existence of conserved moieties in the system, the rate equations can be reduced.

This is the rate equation for all states except MKKK—PP, MKK—PP.

In[17]:=

In Mathematica 5 NDSolve can handle some differential algebraic equations; NDSolve[mapk,{t,0,300}] gives the same result and can be used in earlier Mathematica versions.

In[18]:=

Out[18]=

Here is a plot of MKKK—PP, MKK—PP, and MAPK—PP.

In[19]:=

[Graphics:HTMLFiles/mapk_41.gif]

This gives a plot of MAPK, MAPK—P, MAPK—PP, and their sum.

In[21]:=

Plot[Evaluate[{MAPK[t], MAPK—P[t], MAPK—PP[t], MAPK[t] + MAPK—P[t] + MAPK—PP[t]}/.sol], {t, 0, 300}, PlotStyle→ {Red, Blue, Green, DarkGreen}] ;

[Graphics:HTMLFiles/mapk_43.gif]

An advantage of having the system in reduced form is that its steady state can be computed algebraically. For this example, the rate expressions all are rational functions of the state. Hence, NSolve can be used to find all steady-state solutions.

Here are the equations for the steady state, with functions that are explicitly dependent on time fixed to their values at t=0.

In[22]:=

$sseqn = reqn/.ParameterRules[mapk]/.totRules/.{x_ '[t] →0, t→0} ;$

This finds solutions of the steady-state equations (this may take several seconds). There are very many solutions!

In[23]:=

Out[24]=

This gives the subset of solutions that are both real and positive.

In[25]:=

$sssol = Complement[allsssol, Cases[allsssol, {___, a_→Complex[_, _], ___} | {___, a_→b_ ? Negative, ___}]]$

Out[25]=

The states approach these values close to t=100 in the transient simulation above.

In[26]:=

Out[26]=

In this case FindRoot also finds the steady state starting in the initial conditions for the rate equations (which is not a good approximation to the steady state).

In[27]:=

Out[27]=

1 Kholodenko, B. N., A. Kiyatkin, F. J. Bruggeman, et al. "Untangling the Wires: A Strategy to Trace Functional Interactions in Signaling and Gene Networks." Proceedings of the National Academy of Sciences of the United States of America 99, no. 20 (2002): 12841-12846 . [ online version ]

* Note: If you do not own a copy of Mathematica, you will need to download a trial version in order to view the notebook document.