MATCONT SBML Tutorial 1: Using MATCONT for researching biological models in SBML
Pieter Pareit May 29, 2011
1 Introduction
This tutorial will go through the steps needed to work with a biological model in a SBML-file by using MatCont, a graphical MATLAB software package for the interactive numerical study of dynamical systems [1, 2]. An example paper from The Journal of Cell Biology is taken, Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades [6]. The corresponding SBML-file is located and part of Figure 3 from the article is recreated, see Figure 1. On one of the locations where MatCont detects a bifurcation, a further study is done.
a b
0 20 40 80 80 100
[MAPKKL, (nM) Stationary rates (nM/s)
Figure 1: Figure 3 from [6].
2 Installation
• MatCont
Download location: http://sourceforge.net/projects/matcont/files/
• libSBML
Download location: http://sourceforge.net/projects/sbml/files/libsbml/
Installation instructions:
http://sbml.org/Software/libSBML/docs/cpp-api/libsbml-installation.html
• SBML Toolbox
Download location: http://sourceforge.net/projects/sbml/files/SBMLToolbox/
1
Installation instructions: http: //sbml. org/Sof tware/SBMLToolbox/docs/ManualV3. pdf
Get the latest version of MatCont, the file to download is called matcont<version>.zip. Extract it in a folder matcont<version>, no further installation is needed. MatCont can be started by changing the path in MATLAB to the folder matcont<version> and typing matcont at the MATLAB command line.
The libraries HbSBML and SBML Toolbox need to be installed. There are three options, depending on your situation and platform. (1) Compile and install from source, (2) use an installer corresponding to your platform or (3) unzip a folder with prebuild libraries in matcont<version>/SBML.
The first option, compiling and installing from source, can be used when this is the only way to install the software on your system. See the relevant documentation for your platform in the installation instructions.
The second option, using the installer, is the best option. Get the correct installer from the download location of libSBML and install it. During installation, make sure that you additionally select the MATLAB interface. Next, from the download location of the SBML Toolbox, get the setup program and install it. On some systems the setup program has to be run from the MATLAB command window, this has to be done with a command like cd('c:\Program Files\SBML\SBMLToolbox<version>\toolbox\') ; install.
The third option might be the only option when using a system where you don't have administrator rights. It will only work on 32-bits Windows systems. From the MatCont download location, get the zip-file SBML_libs_win32.zip, this zip-file contains a folder called Win32. Place it in matcont<version>/SBML/. Now you should have a folder matcont<version>/SBML/win32.
3   Paper and SBML-file
Download the pdf of the paper at http://www.jcb.org/cgi/doi/10.1083/jcb.200308060.
Many published biological models are maintained in the BioModels Database at http://www.ebi.ac.uk/ biomodels-main/. Go to that website and enter the name of the main author, Markevich, in the search box. A selection of all models with the author's name in it will be shown. Several models are available, the needed model is BioModels ID BIOMD0000000027, this can be checked in the relevant notes. Choose to download the curated SBML-file. The file BI0MD0000000027.xml will be stored on your computer.
BIOMD0000000027 - Markevich2004_MAPK_orderedMM
Download SBML            | |	Other formats (auto generated)		i	Actions
SBMLL2V1 (curated)				
SBML L2 V2 (auto-generated)m SBML L2 V3 (auto-generated)^ SBML L2 V4 (auto-generated) Publication ID: 14744999			Reference Publication	
		J Cell Biol 2004 Feb;164(3):S53-S. Signaling switches and bistability arising from Markevich Nl. Hoek JB. Kholodenko BN Department of Pathology. Anatomy, and Cell E Philadelphia. PA 19107. USA. [morel		
Hod.,				
Oriainal Model: BIOMD00000000Z7.xml.onam		set#1	bqbioliis	Taxonomy Xenopus laevis
				
Submitter: Nicolas Le Novere		set #2	babiol:isVersionOf Gene Onto-lo-av T	
				
4   Working with the SBML-file 4.1   Importing into MatCont
After starting MatCont with the command matcont in MATLAB, the main MatCont window appears. Select the menu 'Select^Systerrwimport SBML', see Figure 2.
2
Figure 2: Main MatCont window with selection to import an SBML-file.
SBML files have the extension .xml. The file selection dialog that opens will filter on files with that extension. Now go to the location where the file BI0MD0000000027.xml was downloaded and click open as in Figure 3.
D	MatCont	
Select Compute	Window   Type   Options Help	*
System Curve Point Type Curve Type Derivative 5 Duration Status
Select SBML file
3mmmwm
	IQMDOOOOOOOOSO.^ml |	] BIQMDQQ0QQ0Q209.xml   i| BIOM DQQ00Q003 ] BIQMDQuOuuOu225.>;rnl   B michaelrc_and_m
	SHUB [	
	IQMDQ00Q00Q088.x-ml |	3 BIOMDQ0OQ0OQ242.xml   i] tvsonZODl.xml
	ID M D ODD ODD D15 4.X ml j	] BiaMa00O00O0249.xrnl
	IQMDQ00Q00Q169.x-ml |	3 BIOMDQ0OQ0OQ257.xml
	ID M D DDD DDD D17 8 X ml j	] BIOMDDDODDOD258.liml
	IQMDQ00Q00Q186.x-ml j	3 BIOMDuc]0uc]0u27Q.xml
	ID M D ODD ODD D195 .X ml j	] BIOHDDDODDOD315.)iml
		
^ Name: EIOI-.ID00O00O002 .'.xml ?5 of Ivpei | (*. xml)
Open   ] |   Cancel |
Figure 3: Select SBML-file in filechooser.
The biological model in the SBML-file is now read. The species, parameters and reactions from the biological model are converted to coordinates, parameters and differential equations for use by MatCont. The System window, Figure 4, will open.
Most importantly, the importer will generate the system of differential equations
M' = v4 ■ Mp - vx ■ M Mp' — vi ■ M — v2 ■ Mp + v3 ■ Mpp - v4 ■ Mp Mpp' — v2 • Mp — v3 • Mpp
in which
MAPKK *klcat
Kml * (M/Kml + Mp/Km2 + 1)
MAPKK * k2cat Km2 * (M/Kml + Mp/Km2 + 1)
MKP3 * k3cat Km3 * (M/Km5 + Mp/KmA + Mpp/'Km3 + 1)
MKP3 * Meat Kmi * (M/Krah + Mp/Kmi + Mpp/Km3 + 1)'
If desirable, modifications can be done. For now, just press 'OK'.
V2
V3
V4
3
Select   Compute   Window   Type   Options Help
Class
System
Curv'i
Paint Type
Curve Type
Derivatives
Duration
ODE
Markevich2Q04_MAPK_ordered
Name
Coordinates Parameters Time
Derivatives
- numerically
- from window
- symbolically
Markř/Kh20Q4 J-l-.Fk_ordřrřdl.ll.l
iMjMn'.Mr.'f
Iklcat.Kmljkícatj Km2Jk?iatJKm?,k4;at,Km4,Km?JMAPKK,M 1st ord       2nd ord 3rc
4th ord CS
5th ord [S
■ M' = (MKP3*Mp*k4íati:':Km4*(MyKrn5 + Mp/Km4 + Mpp/Km3 + 1)) - (M'MA |Mp' = (M*MAPKK*kleat)/(Kml"(M/Kml + Mp/Km2 + 1)) - (MAPICK*Mp*lc2cai' Mr.-p =-l-li.PKK*l-l|:-*k2íat--Km2ř-l.1-Kml + Mp/Krn2 + 1)) - CMKP3*Mpp"k3
Figure 4: System window with imported system.
4.2   Numerical integration until equilibrium
Now choose the menu 'Type ^Initial Point^Point' so that the Starter and Integrator window will open as in Figure 5.
Initial Point	
E>Z	
M |500	
HP |o	
Mpp |0	
klcat |o.oi	
Kml |50	
k2cat |15	
Km2 |500	
k3cat |o.OS4	
Km3	
k4cat |o.os	
Km4 |ls	
Km5 hg	
MAPKK |50	
MKP3 |100	
Select Cycle	-
Integrator
Integration data
Q®
InitStepsize MaxStepsize Rel, Tolerance Abs. Tolerance F.-=fin>=
Normcontrol
<automatic>
<automatic>
Figure 5: Starter and integrator window with initial values.
The initial values are already filled in as retrieved from the model during the import. From Figure 1 it is clear that for [MAPKK]tot = 50 there are three steady states. Starting from [M] = 500 and [Mp] = [Mpp] = 0 we find the first one by a forward integration.
To obtain this, we leave the default method of integration, ode45, in the Integrator window, but change the interval of integration to 2000. To see visually what is happening, we open the numerics window with the menu 'Window^Numerics' and we open a 2Z?-plot with the menu 'Window^Graphics^2Dplot'. Since we are doing a time integration, it is natural to change the abscissa from 'Coordinates' to 'Time'. We are interested in the behavior of the coordinate Mpp, so change the ordinate to coordinate Mpp. Now on the 2Dplot window, choose the menu 'Layout^Plotting region' and make sure that the abscissa runs from 0 to 2000 and the ordinate from 0 to 100.
4
Using forward integration (menu 'Compute^Forward'), [Mpp] will reach its steady state, close to 50. Selecting 'Compute^Extend' a couple of times, the value of [Mpp] will get stuck in its equilibrium. The setup and result should be as in Figure 6.
■	matcont Q(n)[x]
Select Compute	Window   Type   Options Help
Class	ODE
System	MarKevieh20G4_N1APK_ordered
Curve	P_0(3)
Point Type	P
Curve Type	0
Derivatives	SSSNN
Duration	13.2 sees
Status	ready
	Initial Point
t	1°
M	|500
Mp	1°
M p p	|oo
klcat	|o.oi
Kml	|5 0
k2cat	|l5
Km 2	|500
k3cat	|0.084
Km 3	I22
k4cat	|o.os
Km 4	|ia
Km 5	|7S
MAPKK	|5 0
MKP3	|l 00
	
Select Cycle	
■	2Dplot:l Q@®
File   Edit View	Insert   Tools   Desktop   Window   Help   Layout   Plot ^
	;s|^3.0E«<£-|a|BB|«a
0        200      400 600
1000     1200 1400
Method
Interval InitStepsize MaxStepsize Rel. Tolerance Abs. Tolerance Refine
Normcontrol
Integration data ode45
<automatic>
<automatic>
H	Numeric	EJ(5)@
Window		v
t	2000	
	Coordinate?	
1.1	437.7859436	
Mp	12.84365907	
Mpp	49.37039728	
Figure 6: Setup and result of numerical integration to first equilibrium.
This should also be done for other values of [Mpp] close by. Let's change the initial point Mpp to [Mpp] — 80. Because of Formula 1 in [6],
[Mp] = Mtat - [M] - [Mpp], (1) either [M] or either [Mp] needs to be changed. Set [M] — 420. A forward computation gives Figure 7.
4.3   Continuation of the system
From the equilibrium we can start a continuation. First the initial point has to be selected. To do this choose the menu 'Selects Initial Point' and select the last point of the orbit in the window. The Starter window is automatically updated with the values of the selected point and contains the values of Table 1.
time	0	2000	4000	0	2000	4000
M	500	437.785	437.730	420	437.692	437.729
Mp	0	12.843	12.851	0	12.853	12.852
Mpp	0	49.370	49.418	80	49.453	49.417
Table 1: Numeric values after time integration.
Since we have an equilibrium of the system, select 'Type ^Initial Point^Equilibrium'. The Continuer window will open, and in the Starter windowthe parameter of continuation can be chosen by clicking the corresponding button, take MAPKK. The setup should be as in Figure 8.
5
Figure 7: The numeric integration of the system.
	Initial Point
II	|437.63258
Mp	|l2.353441
Mp p	|49.453933
Jklcat	|o.oi
JKml	|50
Jk2cat	|l5
J)Km2	|500
	|0.034
JKm3	1"
Jk4cat	|0.06
jKm4	|is
JKm5	
(• MAPKK	|50
OMKP3	|ioo
	acobian Data
increment	|le-05
H Contiruer	00®	
Continuation Data		
lnitStep5ize	0.01	
M i n 3te p s i z e	le-5	
SI a;; Stepsi ze	0.1	
Corrector D	ata	
MaxNewtonlters	3	
MaxCorrlters	10	
MaxTestlters	10	
VarTolerance	le-5	
FunTolerance	le-S	
TestTolerance	le-5	
Adapt	3	
Stop Dats		
MaxNumPoints	300	
ClosedCurve	50	
H Numeric	□©a	
Window		
Coordinates		
LI		-
Mp		
Mpp		
Parameters		
klcat		
Ifml		
Figure 8: Updated starter window.
When selecting 'Compute^Forward' to start the continuation, we get the following error: 'No convergence at x0'. The problem is that the system is linearly dependent. By (1) the sum of the concentrations remains constant.
Therefore the system has to be edited to make the equations lineary independent, choose 'Select^Systerrw Edit/Load'. In the Systems selection window, choose the current system, Markevich2004_MAPK_orderedMM.m and select 'Edit'. In the System window, change the equations by replacing the equation for Mp' by Mp — 500 — M — Mpp and remove Mp from the 'Coordinates'. The System window should now look like Figure 9.
6
a	Mat Con t	
Selen Compuie	Window   Type   Options Help	«
Class System Curve Point Type Curve Type Derivatives Duration Status
ODE
Markevich20Q4_MAPK_ordered
Coordinates Parameters
■Markevkh2 004_MAPK_orderedMM_I
|M,Mpp
|klcatJKml,k2cat^m2Jk3catJ<mSJk4catJCm4,Krn5JMAPKKJM
Derivatives
- numerically
- from window
- symbolically
J
J
J J (i
3rd ord
J
4th ord 5th ord S (1
|Mp-500-M-Mpp
|M, = (MKP3*Mp*k4catW(Km4*(M/Km5 + Mp/Km4 + Mpp/Km3 + 1» - CM*MA lMpp, = (MAPKK*Mp*k2cat)/(Krn2*<M/Krnl + Mp/Km2 + 1» - (MKP3*Mpp*k3
Figure 9: Edited system to make it linear independent.
The numeric integration has to be redone. Make sure that the starting coordinates are [M] — 500 and [Mpp] — 0. Select 'Type ^Initial point^Point'. Now do a forward computation with 'Compute^Forward'. Again select as initial point the last point ('Select ^Initial point'). The initial point should be [M] — 437.78813 and [Mpp] = 49.382181. Select 'Type^lnitial Point^Equilibrium' and in the Starter window choose MAPKK as the free parameter. The Starter window should now look as Figure 10.
I.I	|437.7SS13
Up p	|49.382181
Jklcat	|o.oi
JKml	|50
Jk2cat	|l5
J)Km2	|500
Ok3cat	|0.034
JKm3	\22
Jk4cat	|0.06
jKm4	|l8
	I7 3
(• MAPKK	|50
OMKP3	|ioo
	Jacobian Data
increment	|le-05
Monitor Singularities
H Continuer		
Continuation	Data	
InitStepsize	0.01	
M i n Ste p s i z e	le-5	
MaxStepsize	0.1	
Corrector D	ata	
MaxNewtonlters	3	
MaxCorrlters	10	
MaxTestlters	10	
VarTolerance	le-S	
FunTolerance	le-S	
TestTolerance	le-5	
Adapt	3	
Stop Data		
MaxNumPoints	300	
ClosedCurve	50	
H Numeric		
Window		
Coordinate?		
II		I
Mpp		
Parameters		
klcat		
		-
Figure 10: Updated starter window of linear independent system.
Open a 2D plot with 'Window^Graphic^2Dplot' and set the abscissa to the parameter MAPKK e [30, 70] and the ordinate to the state variable (coordinate) Mpp e [0, 500]. Use the menus 'Layout^Variables on axes' and 'Layout^Plotting region' on the 2DPIot window. Choosing 'Compute^Forward' or 'Compute^Backward' and 'Compute^Extend' one or more times, will continue the equilibrium. When MatCont detects a bifurcation, it will pause at that position. To resume the continuation, press 'Resume' as in Figure 11.
Figure 12 shows the complete continuation.
Drawing the curve can be sped up by changing 'Options^Output'. By default this value is 1. This can be
7
■	MatCont Qfnjfx]''
Help i	
Class	ODE
System	Markevich2004_MAPK_ordered
Curve	EP_EP(2)
Point Type	EP
Curve Type	EP
Derivatives	SSSNN
Duration	
Status	Limit point
M	437.7SS13
Mpp	49.382181
!>lcat	O.Ol
JKml	50
Jk2<at	15
JKm2	500
Jk3cat	0.084
JKm3	22
Jk4cat	0.06
■ . ,- i	18
JKm5	73
(• MAPKK	50
JMKP3	100
J	acobian Data
increment	|1»-0S
File   Edit   View   Jnsen   Tools   Desktop   Window   Help   Layout Plot
xNumPoints
b	Numeric q@®	
Window		
	Coordinates	
ri	352.21592	
Mpp	129.14652 Parameters	
fclcat	0.01 5n	
Figure 11: Setup and output of an equilibrium continuation.
*
Figure 12: Continuation of equilibrium.
set to a higher value, in Figure 13 this is set to 100. While the calculations will be as precise as before, the drawing will now only be done after 100 steps are calculated.
100 P°in,!
°K I r-™«<
Figure 13: Changing the interval between redraws.
By default, the continuer will continue 300 points. In this case more than 5000 will have to be continued. After the first 300 points, we can extend the curve with 'Compute^Extend'. If it is known that the curve needs to be extended many times, MaxNumPoints in Continuer window can be set to a larger value.
8
4.4   Continuation in 2 free parameters
Continuation in 2 parameters is possible. Then the bifurcation boundaries corresponding to the limit point bifurcations are curves in the parameter plane. For instance, select the second limit point with 'Select ^Initial point'. MatCont will now do a LP-LP continuation. The starter loads automatically the values in and will look as in Figure 14. We select kml and MAPKK as free parameters.
	Initial Point
II	|4S. 763794
Mpp	|43S.452S3
BP Active	
J Jklcat	|o.oi
J (iKml	1 I50
J Jk2cat	|l5
J JKm2	|500
J Jk3cat	|0.0S4
J jKm3	|z2
J Jk4cat	|0.06
J JKml	|l8
J jKm5	l7S
J   (i MAPKK	|39.246024
J JMKP3	100
Jacobian Data	
increment	le-05
W             Conti nuer	Q@©1	
Continuation Data		
lnitStepsiz:e	0.01	
MinStepsi&e	le-5	
Ma;:^tep>si;ze	0.1	
Corrector Data		
MaxNewtonlters	3	
MaxCorrlters	10	
MaxTestlters	10	
VarTolerance	U-6	
FunTolerance	le-6	
TestTolerance	le-5	
Adapt	3	
Stop DatE		
MaxNumPoints	5000	
CloseclCurva	50	
■	Numeric [T][H][xf	
Window		
	Coordinate;	
|.|	46.9625	
Mpp	451.2744	
	Parameter;	
klcat	0.00019207673	
Kml	5fi	
Figure 14: Starter window for LP-LP continuation.
It is best to look at the parameter plane. Select a 2DPIot with on the axes the relevant parameters, MAPPK and Kml. A good plotting region for MAPKK is from 35 to 60 and for Kml from 30 to 270. A forward continuation ('Compute^Forward') gives Figure 15. At M = 249.921974, Mpp = 240.426057, kml = 246.269558, MAPKK = 41.157105 a cusp bifurcation is detected.
Figure 15: LP-LP continuation.
9
References
[1] A. Dhooge, W. Govaerts, and Yu.A. Kuznetsov. MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs. ACM Transactions on Mathematical Software (TOMS), 29(2):141-164, 2003.
[2] A. Dhooge, W. Govaerts, Yu.A. Kuznetsov, W. Mestrom, A.M. Riet, and B. Sautois. MATCONT and CL_MATCONT: Continuation toolboxes in Matlab. http://sourceforge.net/projects/matcont/, 2004.
[3] W. Govaerts, Yu.A. Kuznetsov, and B. Sautois. MATCONT. Scholarpedia, 1(9):1375, 2006.
[4] M. Hucka, A. Finney, H.M. Sauro, H. Bolouri, J.C. Doyle, H. Kitano, A.P. Arkin, B.J. Bornstein, D. Bray, A. Cornish-Bowden, A.A. Cuellar, S. Dronov, E.D. Gilles, M. Ginkel, V. Gor, I.I. Goryanin, W.J. Hedley, T.C. Hodgman, J.H. Hofmeyr, P.J. Hunter, N.S. Juty, J.L. Kasberger, A. Kremling, U. Kummer, N. Le Novere, L.M. Loew, D. Lucio, P. Mendes, E. Minch, E.D. Mjolsness, Y. Nakayama, M.R. Nelson, P.F. Nielsen, T. Sakurada, J.C. Schaff, B.E. Shapiro, Shimizu T.S., H.D. Spence, J. Stelling, K. Takahashi, M. Tomita, J. Wagner, and J. Wang. The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics, 19(4):524-531, 2003.
[5] C. Li, M. Donizelli, N. Rodriguez, H. Dharuri, L. Endler, V. Chelliah, L. Li, E. He, A. Henry, M.I. Stefan, etal. Biomodels database: An enhanced, curated and annotated resource for published quantitative kinetic models. BMC systems biology, 4(1):92, 2010.
[6] N.I. Markevich, J.B. Hoek, and B.N. Kholodenko. Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades. The Journal of cell biology, 164(3):353, 2004.
10