This is a continuation of last week's post. I had some troubles getting the graph that I wanted, but Bard Ermentrout (the creator of XPP) fixed the problem in a sentence via email (lesson: when stuck, don't be afraid to ask the experts). There will be many technical terms which won't be elaborated on in this post. I encourage the reader who want to follow along to look up things they are not familiar with on Wikipedia. The post assumes you have some familiarity with bifurcation analysis (the qualitative study of ODEs).
Code 1.
A' = - kas1*A*I + kdi1*C + kcat1*C
IM' = - km*A*IM - kdm1*IM + kdm2*IMM + kcat1*C
IMM'= km*A*IM - kdm2*IMM
C = AT-A
I = IT-C-IM-IMM
par IT=0.01, AT=1, kas1=100, kdi1=0.5, km=50
par kcat1=0.5, kdm1=1, kdm2=1
@ total=500,dt=0.25,meth=STIFF,xlo=0,xhi=500,ylo=0,yhi=2
@ NTST=12,NMAX=100000,NPR=100000,DS=0.01
@ DSMAX=0.01,DSMIN=0.0001,PARMIN=0,PARMAX=10
@ AUTOXMIN=0,AUTOXMAX=10,AUTOYMIN=0,AUTOYMAX=4
done
Code 2.
AI = 2*AT*IH/(AT+IH+Km1 +sqrt((AT+IH+Km1)^2-4*AT*IH))
A = AT - AI
IM = IT - IH - IMM
dIH/dt = kdm1*IM - kcat1*AI
dIMM/dt = kas2*A*IM - kdm2*IMM
par AT=2, IT=6, Km1=0.01
par kcat1=0.5, kas2=50, kdm1=1, kdm2=1
# use AT=3 for Fig.3D
@ xp=IH, yp=IMM
@ total=400,dt=0.5,meth=STIFF
@ xlo=0,xhi=6,ylo=0,yhi=6,NMESH=400
done
This code comes from the supplement of a simulation made by Verdugo et al. last year [1]. The context is modeling cell cycle checkpoints as bistable switches. It is a very exciting model and framework, building on decades of work [2], but because of time-constrains we won't get into the details of it in this post.
Load code 1. Start by running the simulation (press I then G), and finding a steady steady state (S G). We find a steady state at around A=1. Then we go to the AUTO software and run for the steady state (R S). We can the grab points of interests (G and tab). For example we see that the second and third points are of type LP, which is stands for limit points and is a saddle-node bifurcation. The trajectories that are red are stable, others not.
A saddle-node bifurcation point is basically a point where you from monostable (one stable fixed point) to bistabile (two stable points and one unstable "decision point" fixed point). We can see this on the nullclines (N N) of the phase plane. We will do this using code 2 later (it is the same basic model, except for some technical simplifications and different parameters).
We grab the second point (press enter once you've selected a point) and then we want to perform a two-parameter bifurcation. Change to two parameter by Axes - Two Param (A T). We have selected IT as our main parameter and AT as our secondary parameter. We run it again and should get a small line that represents that bifurcation point in the previous picture, but now also depending on the second parameter. So for any given AT, like for example 1, we see two bifurcation points, which correspond to the previous diagram.
We have to repeat this process a couple of times - first for the other limit point 3, then to go "backwards" in the other direction for numerical simulation reasons (thanks again to Bard Ermentrout for pointing this out to me - I'm sure this is obvious to people who have done more numerical analysis, but it wasn't obvious to me). This is done by first selecting run and go to the last 1-par graph (Axes menu) to get the 1-parameter bifurcation (press D for redraw if it doesn't show up), then select the new point, then run with the last 2-par diagram. Then change the Ds, "delta", under the Numerics tab to go in the other direction, i.e. add a minus sign in front. Once both directions for both points have been covered it should look like this.
The way to interpret this is that the region in the middle is the bistabile region, that is where cell cycle checkpoints can occur. The region above that is where the checkpoint is disengaged, and the one below where it is engaged. If you look at the first picture carefully, and maybe by changing the parameter AT and seeing how the shape changes, you can convince yourself that this is true. For example, a low IT, say 2, and AT=1 will correspond to a high A in the first picture, which corresponds to us being in the upper left region and thus a disengaged checkpoint - all systems go.
Here are the null clines produced by code 2. In the first picture AT is 2 and in the second 3. We can here see the bifurcation point in action - in the first case it is bistabile, in the second monostabile. The circles correspond to stabile fixed points and the triangles to unstable fixed points. A simple way to get them automagically is to use a basic Monte Carlo simulation (S C).
That's All Folks!
1: Here's the paper http://rsob.royalsocietypublishing.org/content/3/3/120179.short
2: For example via the Tyson Lab http://mpf.biol.vt.edu/lab_website/