Bifurcation diagram How-To
745 views
Skip to first unread message
Pedro Mendes
May 27, 2015, 6:25:20 PM5/27/15
to copasi-u...@googlegroups.com
THIS IS A LONG POST. TLDR: it goes through all the steps to create a
bifurcation diagram. Eventually a video will be made following this
description, but for now I'm afraid that text was all that was possible!
HOW TO DO BIFURCATION ANALYSIS WITH COPASI
In bifurcation analysis we want to observe how the behaviour of the
system changes while changing one parameter -- named the bifurcation
parameter. In this type of nonlinear analysis we want to see how the
system changes from one stable steady state, to an oscillation, or
between oscillations with different periodicity, or even transitions to
higher complexity attractors (such as chaos). If the system only has
steady states (no oscillations) then this type of analysis is useless...
Here are the steps that you should take to do a bifurcation analysis of
a model.
1- load your model. In this example we are going to use the Kholodenko
2000 model of MAPK cascade. You can find this in the curated branch of
the BioModels database:
https://www.ebi.ac.uk/biomodels-main/BIOMD0000000010
Download the SBML file and import it into COPASI.
2- In COPASI bifurcation diagrams are constructed by iterating
cross-sections. Cross-sections are hyper-surfaces in phase-space that
are crossed by the time course trajectory. These hypersurfaces are
selected by picking one of the variables of the model (usually a
species) and deciding on a specific value of that variable. Everytime
that the trajectory passes through that value of the variable, in one
predetermined direction (negative or positive), a point is plotted on
that hypersurface.
3- So in the model of interest we need pick a) a variable, b) a value
for that variable and c) a direction. It is important that the value be
one that the variable takes during a trajectory. The best way to chose
this is to first run a simple time course and observe the variable of
interest. Chose a value taken by the variable (it could be any value
that it takes during the time course). In the MAPK model let's chose the
MEK1-PP species and the value 10 nM; for direction we pick the positive
(ie when the values are increasing).
4- To visualize the cross-section we create a 2D phase-space plot, ie a
plot that has a variable on the X axis and another on the Y-axis. Let's
create it for the pair of species ERK2-PP and MOS-P (follow the video
tutorial on creating plots if you don't know how to do this). It is
important that you pick the transient concentrations of the species and
that any of the species you chose for the plot not be the one you chose
above for the cross section (ie not MEK1-PP). Make sure that you select
Symbols for the plot (could also be points, but not lines).
5- Since cross-sections are analyses of time courses, let's first run a
time course. Go to the time course page, select a long-enough end time
(ie you want to see the whole behaviour, in this case several
oscillations) -- let's use use 30000, and also set 600 intervals-- and
press run. You should see a sequence of points in the phase plane,
orienting along a circular line (the limit cycle). Good, now we're ready
for the cross section
6- Go to the "Cross Section" task and chose the chosen variable for the
cross section (transient concentration of MEK1-PP), set the "crossing"
value as 10, and chose positive direction. We have to tell the task that
it can collect crossings until the end time of the time course (in this
case 30000). This number needs to be entered in the box labelled "Stop
detecting crossings after any condition" "if detection time(s) larger".
Now press run and observe the phase space diagram.
7- The phase space diagram shows two points. Each point corresponds to
an oscillation that passed the value of ERK1-PP=10nM when the variable
was increasing. Actually, if you look closesly, there are more than two
points, but they are too close to be distinguished. If you zoom in
around the point on the lower right you will see that there are actually
two points here. Because this system oscillates, there are several
crossings. The first one corresponds to the first oscillation that was
larger than the others. The remaining ones are all very close (you saw
two of them when you zoomed in), in fact many are overlayed on top of
each other. Since they are very close, this means that you have a stable
oscillation. If you had a simple steady state, then the phase plane
would be empty as in a steady state no concentrations change and so
ERK1-PP does not "pass through" 10nM during the steady state. In
cross-section analysis steady states are invisible. A simple (period-1)
oscillation results in one point, a complex oscillation results in
several points (eg period-2 by two points, period-3 by 3, etc) and chaos
will be a cloud of many points.
8- To clean up this cross-section plot, and obtain only one point
corresponding to the stable oscillation, we need to get rid of the first
few oscillations that occurred when the system was still stabilizing.
Let's do this by not collecting any crossings until time=9000. For this
type 9000 in the box that says "Find crossings after any condition is
met" "if time (s) is larger than". This means that the crossings are
only going to be plotted after 9000 second has passed and this should be
enough to remove the early oscillations. Run again. Now you should see
only one point in the phase-space diagram because this is a simple
period-1 oscillation. Actually there were more than one crossings, but
they happened all at the same values of MOS-P and ERK2-PP, so we only
see one point.
Now that we finally created a cross-section diagram for this oscillatory
system, we are ready to create a bifurcation diagram. This will be a
diagram that shows how the oscillations change when we change a
parameter: the bifurcation parameter.
9- So now we need to find such a bifurcation parameter. In this model
let's use the parameter Ki in the reaction "MAPKKK activation" (this
parameter represents the strength of feedback inhibition of this
reaction by the ERK2-PP species).
10- We need to create a new plot to visualize the bifurcation diagram.
This plot should have our parameter (Ki of reaction "MAPKKK activation")
in the X axis, and some variable (not the one chosen for the
cross-section) in the Y axis. Let's pick MOS-P. Create the new plot and
call it "bifurcation diagram"). I create the plot using points, as we
are going to have many iterations.
11- The bifurcation diagram is essentially a parameter scan of a
cross-section. So we now go to the "Parameter Scan" task and create a
parameter scan for Ki of reaction "MAPKKK activation" (it is in
reactions-"MAPKKK activation"-kinetic parameters-Ki). Set the task to be
"Cross-Section", make sure that "output during subtask execution" is on.
Make the scan for 500 intervals between 0.01 and 50 (you may want it
logarithmic, but it is ok to have it non-log as well). Run and observe
the bifurcation diagram.
12- A few important things come out of the bifurcation diagram. Each
point in the diagram corresponds to a different stable oscillation.
These oscillations have different properties (the concentration of MOS-P
when the oscillation passes ERK-PP=10nM) but they change rather smoothly
with the bifurcation parameter (on the X axis). Interestingly, even
though the parameter varied all the way to 50, there are no more points
for Ki>20. If you change the plot to log scale (on the menu of the plot
window, press "Log X"), you will see that even though the parameter
started at Ki=0.01, the first point appears at about Ki>0.1. This means
that there are only stable oscillations inside the range 0.1<Ki<20. At
Ki~=0.1 there is a bifurcation when the system goes from a steady state
to stable (period-1) oscillations. At Ki~=20 there is another
bifurcation when the system transitions to a steady state.
13- To check this, set the value of the parameter Ki of reaction "MAPKKK
activation" in the model to 25 (ie beyond the second bifurcation). Then
go to the Time Course task and run (you should create a time course plot
if you haven't done yet). You will see that the dynamics now is of a
steady state, with some damped oscillations in the transient. Now change
the value to 0.01 and run again. You see a similar thing, except in this
case the damped oscillations have higher frequency. You can also create
a slider for this parameter in the Time Course task, set its boundaries
to be between 0.01 and 50, and manually change it and see how the
dynamics changes. Around 20 you will see the oscillations becoming
damped rather than stable -- this is a Hopf bifurcation.
The bifurcation diagram is a convenient way to check the domain of
parameter values in which there are oscillations or more complex
behaviour. This example model has only simple, period-1, oscillations
but in other cases you may observe the line bifurcating into two lines
(a bifurcation changing the oscillation from period-1 to period-2), and
maybe those two into several others (what is sometimes called "the route
to chaos").
There are other ways to construct bifurcation diagrams, but currently
this is the only way in which you can do it in COPASI: scanning the
bifurcation parameter and plotting cross sections at each value.
Happy bifurcations!
--
Pedro Mendes
Professor of Computational Systems Biology
School of Computer Science
Manchester Centre for Integrative Systems Biology
University of Manchester
Manchester Institute of Biotechnology
131 Princess Street
Manchester, M1 7DN, U.K.
bifurcation diagram. Eventually a video will be made following this
description, but for now I'm afraid that text was all that was possible!
HOW TO DO BIFURCATION ANALYSIS WITH COPASI
In bifurcation analysis we want to observe how the behaviour of the
system changes while changing one parameter -- named the bifurcation
parameter. In this type of nonlinear analysis we want to see how the
system changes from one stable steady state, to an oscillation, or
between oscillations with different periodicity, or even transitions to
higher complexity attractors (such as chaos). If the system only has
steady states (no oscillations) then this type of analysis is useless...
Here are the steps that you should take to do a bifurcation analysis of
a model.
1- load your model. In this example we are going to use the Kholodenko
2000 model of MAPK cascade. You can find this in the curated branch of
the BioModels database:
https://www.ebi.ac.uk/biomodels-main/BIOMD0000000010
Download the SBML file and import it into COPASI.
2- In COPASI bifurcation diagrams are constructed by iterating
cross-sections. Cross-sections are hyper-surfaces in phase-space that
are crossed by the time course trajectory. These hypersurfaces are
selected by picking one of the variables of the model (usually a
species) and deciding on a specific value of that variable. Everytime
that the trajectory passes through that value of the variable, in one
predetermined direction (negative or positive), a point is plotted on
that hypersurface.
3- So in the model of interest we need pick a) a variable, b) a value
for that variable and c) a direction. It is important that the value be
one that the variable takes during a trajectory. The best way to chose
this is to first run a simple time course and observe the variable of
interest. Chose a value taken by the variable (it could be any value
that it takes during the time course). In the MAPK model let's chose the
MEK1-PP species and the value 10 nM; for direction we pick the positive
(ie when the values are increasing).
4- To visualize the cross-section we create a 2D phase-space plot, ie a
plot that has a variable on the X axis and another on the Y-axis. Let's
create it for the pair of species ERK2-PP and MOS-P (follow the video
tutorial on creating plots if you don't know how to do this). It is
important that you pick the transient concentrations of the species and
that any of the species you chose for the plot not be the one you chose
above for the cross section (ie not MEK1-PP). Make sure that you select
Symbols for the plot (could also be points, but not lines).
5- Since cross-sections are analyses of time courses, let's first run a
time course. Go to the time course page, select a long-enough end time
(ie you want to see the whole behaviour, in this case several
oscillations) -- let's use use 30000, and also set 600 intervals-- and
press run. You should see a sequence of points in the phase plane,
orienting along a circular line (the limit cycle). Good, now we're ready
for the cross section
6- Go to the "Cross Section" task and chose the chosen variable for the
cross section (transient concentration of MEK1-PP), set the "crossing"
value as 10, and chose positive direction. We have to tell the task that
it can collect crossings until the end time of the time course (in this
case 30000). This number needs to be entered in the box labelled "Stop
detecting crossings after any condition" "if detection time(s) larger".
Now press run and observe the phase space diagram.
7- The phase space diagram shows two points. Each point corresponds to
an oscillation that passed the value of ERK1-PP=10nM when the variable
was increasing. Actually, if you look closesly, there are more than two
points, but they are too close to be distinguished. If you zoom in
around the point on the lower right you will see that there are actually
two points here. Because this system oscillates, there are several
crossings. The first one corresponds to the first oscillation that was
larger than the others. The remaining ones are all very close (you saw
two of them when you zoomed in), in fact many are overlayed on top of
each other. Since they are very close, this means that you have a stable
oscillation. If you had a simple steady state, then the phase plane
would be empty as in a steady state no concentrations change and so
ERK1-PP does not "pass through" 10nM during the steady state. In
cross-section analysis steady states are invisible. A simple (period-1)
oscillation results in one point, a complex oscillation results in
several points (eg period-2 by two points, period-3 by 3, etc) and chaos
will be a cloud of many points.
8- To clean up this cross-section plot, and obtain only one point
corresponding to the stable oscillation, we need to get rid of the first
few oscillations that occurred when the system was still stabilizing.
Let's do this by not collecting any crossings until time=9000. For this
type 9000 in the box that says "Find crossings after any condition is
met" "if time (s) is larger than". This means that the crossings are
only going to be plotted after 9000 second has passed and this should be
enough to remove the early oscillations. Run again. Now you should see
only one point in the phase-space diagram because this is a simple
period-1 oscillation. Actually there were more than one crossings, but
they happened all at the same values of MOS-P and ERK2-PP, so we only
see one point.
Now that we finally created a cross-section diagram for this oscillatory
system, we are ready to create a bifurcation diagram. This will be a
diagram that shows how the oscillations change when we change a
parameter: the bifurcation parameter.
9- So now we need to find such a bifurcation parameter. In this model
let's use the parameter Ki in the reaction "MAPKKK activation" (this
parameter represents the strength of feedback inhibition of this
reaction by the ERK2-PP species).
10- We need to create a new plot to visualize the bifurcation diagram.
This plot should have our parameter (Ki of reaction "MAPKKK activation")
in the X axis, and some variable (not the one chosen for the
cross-section) in the Y axis. Let's pick MOS-P. Create the new plot and
call it "bifurcation diagram"). I create the plot using points, as we
are going to have many iterations.
11- The bifurcation diagram is essentially a parameter scan of a
cross-section. So we now go to the "Parameter Scan" task and create a
parameter scan for Ki of reaction "MAPKKK activation" (it is in
reactions-"MAPKKK activation"-kinetic parameters-Ki). Set the task to be
"Cross-Section", make sure that "output during subtask execution" is on.
Make the scan for 500 intervals between 0.01 and 50 (you may want it
logarithmic, but it is ok to have it non-log as well). Run and observe
the bifurcation diagram.
12- A few important things come out of the bifurcation diagram. Each
point in the diagram corresponds to a different stable oscillation.
These oscillations have different properties (the concentration of MOS-P
when the oscillation passes ERK-PP=10nM) but they change rather smoothly
with the bifurcation parameter (on the X axis). Interestingly, even
though the parameter varied all the way to 50, there are no more points
for Ki>20. If you change the plot to log scale (on the menu of the plot
window, press "Log X"), you will see that even though the parameter
started at Ki=0.01, the first point appears at about Ki>0.1. This means
that there are only stable oscillations inside the range 0.1<Ki<20. At
Ki~=0.1 there is a bifurcation when the system goes from a steady state
to stable (period-1) oscillations. At Ki~=20 there is another
bifurcation when the system transitions to a steady state.
13- To check this, set the value of the parameter Ki of reaction "MAPKKK
activation" in the model to 25 (ie beyond the second bifurcation). Then
go to the Time Course task and run (you should create a time course plot
if you haven't done yet). You will see that the dynamics now is of a
steady state, with some damped oscillations in the transient. Now change
the value to 0.01 and run again. You see a similar thing, except in this
case the damped oscillations have higher frequency. You can also create
a slider for this parameter in the Time Course task, set its boundaries
to be between 0.01 and 50, and manually change it and see how the
dynamics changes. Around 20 you will see the oscillations becoming
damped rather than stable -- this is a Hopf bifurcation.
The bifurcation diagram is a convenient way to check the domain of
parameter values in which there are oscillations or more complex
behaviour. This example model has only simple, period-1, oscillations
but in other cases you may observe the line bifurcating into two lines
(a bifurcation changing the oscillation from period-1 to period-2), and
maybe those two into several others (what is sometimes called "the route
to chaos").
There are other ways to construct bifurcation diagrams, but currently
this is the only way in which you can do it in COPASI: scanning the
bifurcation parameter and plotting cross sections at each value.
Happy bifurcations!
--
Pedro Mendes
Professor of Computational Systems Biology
School of Computer Science
Manchester Centre for Integrative Systems Biology
University of Manchester
Manchester Institute of Biotechnology
131 Princess Street
Manchester, M1 7DN, U.K.
Pnar Pir
May 27, 2015, 6:47:56 PM5/27/15
to copasi-u...@googlegroups.com
Thanks Pedro, it will treasure this email for future use...
Recently I was pointed to this paper for identifiability analysis in COPASI and it saved me so much time: http://www.med.uni-magdeburg.de/fme/sysbio/wp-content/uploads/papers/BIO3329.pdf
I think publishing such short 'method' papers is a great idea... Maybe Copasi.org can catalogue such papers in addition to video tutorials?
Best,
Pnar
--
You received this message because you are subscribed to the Google Groups "COPASI User Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to copasi-user-fo...@googlegroups.com.
To post to this group, send email to copasi-u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/copasi-user-forum/55664477.7070406%40manchester.ac.uk.
For more options, visit https://groups.google.com/d/optout.
Pnar Pir
----
A mistake repeated more than once is a decision - Paulo Coelho
Donrad
Dec 15, 2016, 12:16:47 AM12/15/16
to COPASI User Forum
Hi All,
First of all thanks for this comprehensive tutorial. I tried it step by step and I have a confusing issue with the output results.For the first time that I run the bifurcation analysis I see the graph that matches to the explained example here. but when I try to run it afterwards (click the run button again i.e. after the first run) the bifurcation plot turns into a series of points oriented all vertically (i.e. toward the bifurcation parameter value). But when I hit run button in the "Time course" task section it shows the very first shape of the graph when i first run the analysis according to the explanation here. I'm confused what's going on?
I'll appreciate if you let me know which one is the correct approach
Reply all
Reply to author
Forward
0 new messages