Computer Laboratory on
Bursting Oscillations in Pancreatic Beta Cells
Workshop on Mathematical Physiology
June 14 - 25, 1999
University of British Columbia
We will be exploring a simple model of bursting electrical activity.
The model is biophysically based, and similar in concept to the Chay-Keizer
model of bursting that was discussed in the introductory lecture on bursting.
The model incorporates the minimal ingredients needed to obtain a burster,
The first two components interact to produce the action potentials during
the active phase, and the last two components interact to switch the system
between the active and silent phases. Comparing this model to the
Chay-Keizer model, s is analogous to free intracellular calcium concentration
- s increases during the active phase, and decreases during the silent
ICa: a fast, voltage-gated, inward current
IK: a slower voltage-gated outward current
Is: an s-sensitive inward current (activates with increasing
s: a slow variable
A) Getting used to the model:
B) Making a bifurcation diagram for the fast subsystem of the model:
Download the model (bmb.ode)
to your workspace, and take a close look at the file. Can you determine
the purpose of the autoflag and sknot parameters?
Start XPP, and compute a numerical solution with the model as is.
In the main window, or in a new window (use the Makewindow menu command),
plot the solution in the (s,v) phase plane (use the Viewaxes menu command).
Also plot the solution in the three-dimensional phase space. Are
the graphs as expected?
We will use the slow variable s (which is equal to sknot under steady
state conditions) as the bifurcation parameter.
C) Creating (v,n) phase planes at strategically chosen values of s:
It is easiest to start the bifurcation diagram at a steady state.
Set the autoflag parameter to 1, and compute another numerical solution.
You should obtain a high-voltage steady state solution, corresponding to
the fixed point of the fast subsystem at s=sknot. Use the Integrate-Last
commands a few times to ensure that you are truly at the steady state.
Start the AUTO component of XPP (use the File-Auto menu command).
You should see a new window appear - this is where the bifurcation diagram
will be plotted. Note that the axes labels are already correctly
set - by default, AUTO uses the first parameter listed in the bmb.ode file
as the bifurcation parameter (note that you can choose a different bifurcation
parameter with the Parameter menu command). To set up the axes, choose
Axes-HiLo, and set Xmin = -0.5, Ymin = -80, Xmax = 1, and Ymax = 0.
To set up the numerics, choose Numerics, and set Nmax = 1000, NPr = 100,
Dsmin = 0.000001, Dsmax = 0.1, ParMin = -2, and ParMax = 2.
Choose Run-SteadyState, and watch the Z-curve appear (a heavy curve indicates
a solid steady state, and a thin curve indicates an unstable steady state).
Note the printout of special points in the window from which you started
XPP (EP stands for End Point; HB for Hopf Bifurcation, LP for Limit Point,
and here indicates a Saddle-Node bifurcation, or fold).
Choose Grab, and use the Tab and arrow keys on your keyboard to navigate
along the bifurcation curve. Note the changing printout at the bottom
of the AUTO window. To get out of Grab-mode, press enter.
Grab the Hopf bifurcation point, choose Run-Periodic, and watch the branch
of periodics appear. AUTO will have trouble near the homoclinic bifurcation
(why?). MX indicates no convergence.
Choose Grab again, and navigate along the branch of periodics. Take
special note of the value of the period (printed at the bottom right-hand
corner of the AUTO window) as the homoclinic bifurcation is approached.
Is it as expected?
Record the values of s at the Hopf bifurcation (sHB), both saddle-nodes
(sSNleft, and sSNright), and the homoclinic bifurcation
(sHC - estimate this value). You will use these values
Save the diagram (choose File-WritePts). You will use the diagram
At each value of s, there is a corresponding phase plane for the fast
subsystem. The phase plane changes qualitatively at the values of
s recorded in B), sHB ,sSNleft ,sHC ,
D) Superimposing the bifurcation diagram and the numerical solution
to the full burster:
Return to the main XPP window. Use the ViewAxes-2D command to set
up the (n,v) phase plane (Xmin = -0.1, Ymin = -80, Xmax = 0.5, Ymax = -10).
Change the total integration time from 14000 to 1000 (use the nUmerics-Total
menu command), and Dt from 5 to 0.1 (use the nUmerics-Dt menu command).
Set s = sknot < sHB, and use several Integrate-Mouse commands
to compute a phase portrait for this value of s. Once you have an
idea about the location of a fixed point, you can use the SingPts menu
command to determine its exact location. Are the results as expected?
Are they consistent with the information obtained from the bifurcation
Make a phase portrait for sHB < s < sSNleft,
and determine the location and stability of any fixed points. Again,
are the results consistent with the bifurcation diagram?
Repeat for sSNleft < s < sHC, sHC
< s < sSNright, and s > sSNright.
E) Investigating the effect of taus:
Set the autoflag parameter back to 0, Total back to 14000, and Dt back
to 5. Use the ViewAxes-2D command to set up the (s,v) phase plane
(Xmin = 0.15, Ymin = -70, Xmax = 0.2, Ymax = -10). Compute a solution
of the full burster (you may have to erase the solution and reintegrate
once or twice to eliminate a transient solution).
To plot the bifurcation diagram you saved in B), choose GraphicStuff-Freeze-BifDiag.
The parameter taus determines the time scale for the variable s.
Since taus is large, s is slow.
F) Investigating the effect of R:
Change the total time of integration to 30000, and change the value of
the parameter taus from 35000 to 75000. To fully observe the effect
of changing taus, you should use two display windows: one as described
under D), and another one that will graph v versus t (Xmin = 0, Ymin =
-70, Xmax = 30000, Ymax = -10). Calculate a numerical solution.
Repeat with taus = 200000. Can you explain what is happening?
Set taus = 10000, and compute another solution. Repeat with taus
= 5000. What is happening and why?
Biophysically, the best candidate for an electrical glucose sensor is
the ATP-sensitive K+ channel (one of the products of glucose
metabolism is ATP, and ATP shuts off the K(ATP) channels, thus depolarizing
the membrane potential). Included in the model is a K(ATP) current,
ikatp, dependent on the parameter R. That is, varying R is similar
to varying the glucose concentration.
Change R from 0.5 to 0.7, and compute a numerical solution. Observe
the solution in the window displaying v versus t. How does the solution
differ from before? Repeat with R = 0.8 and R = 0.9. Can you
explain what is happening (in terms of the K(ATP) current)?
Repeat with values of 0 < R < 0.5. Again, what is happening,
and can you explain why?
Challenge: How does the bifurcation diagram depend on R? Can
you explain the effect of R on the electrical activity in terms of changes
in the bifurcation diagram?