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 ChayKeizer
model of bursting that was discussed in the introductory lecture on bursting.
The model incorporates the minimal ingredients needed to obtain a burster,
namely

I_{Ca}: a fast, voltagegated, inward current

I_{K}: a slower voltagegated outward current

I_{s}: an ssensitive inward current (activates with increasing
s)

s: a slow variable
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
ChayKeizer model, s is analogous to free intracellular calcium concentration
 s increases during the active phase, and decreases during the silent
phase.
A) Getting used to 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 threedimensional phase space. Are
the graphs as expected?
B) Making a bifurcation diagram for the fast subsystem of the model:
We will use the slow variable s (which is equal to sknot under steady
state conditions) as the bifurcation parameter.

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 highvoltage steady state solution, corresponding to
the fixed point of the fast subsystem at s=sknot. Use the IntegrateLast
commands a few times to ensure that you are truly at the steady state.

Start the AUTO component of XPP (use the FileAuto 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
AxesHiLo, 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 RunSteadyState, and watch the Zcurve 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 SaddleNode 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 Grabmode, press enter.

Grab the Hopf bifurcation point, choose RunPeriodic, 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 righthand
corner of the AUTO window) as the homoclinic bifurcation is approached.
Is it as expected?

Record the values of s at the Hopf bifurcation (s_{HB}), both saddlenodes
(s_{SNleft}, and s_{SNright}), and the homoclinic bifurcation
(s_{HC } estimate this value). You will use these values
in C).

Save the diagram (choose FileWritePts). You will use the diagram
in D).
C) Creating (v,n) phase planes at strategically chosen values of s:
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), s_{HB} ,s_{SNleft} ,s_{HC} ,
and s_{SNright}.

Return to the main XPP window. Use the ViewAxes2D 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 nUmericsTotal
menu command), and Dt from 5 to 0.1 (use the nUmericsDt menu command).
Set s = sknot < s_{HB,} and use several IntegrateMouse 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
diagram?

Make a phase portrait for s_{HB} < s < s_{SNleft},
and determine the location and stability of any fixed points. Again,
are the results consistent with the bifurcation diagram?

Repeat for s_{SNleft} < s < s_{HC}, s_{HC}
< s < sS_{Nright}, and s > s_{SNright}.
D) Superimposing the bifurcation diagram and the numerical solution
to the full burster:

Set the autoflag parameter back to 0, Total back to 14000, and Dt back
to 5. Use the ViewAxes2D 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 GraphicStuffFreezeBifDiag.
Cool, eh.
E) Investigating the effect of taus:
The parameter taus determines the time scale for the variable s.
Since taus is large, s is slow.

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?
F) Investigating the effect of R:
Biophysically, the best candidate for an electrical glucose sensor is
the ATPsensitive 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?