Example: A MAPK Cascade
This model is taken from B. N. Kholodenko et al. 2002 ^{1}. 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.
A mitogenactivated 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 righthand side of the mass balance equations.
In[7]:=
Out[7]=
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]:=
Out[15]=
This gives a list of rules for the total amounts of MKKK, MKK, and MAPK.
In[16]:=
Out[16]=
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]:=
This gives a plot of MAPK, MAPK—P, MAPK—PP, and their sum.
In[21]:=
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 steadystate 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]:=
This finds solutions of the steadystate 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]:=
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):
1284112846 . [ 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.
