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, namely

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 phase.

A) Getting used to the model:

  1. 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?
  2. 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?
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.

  1. 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.
  2. 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.
  3. 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).
  4. 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.
  5. 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.
  6. 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?
  7. 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 in C).
  8. Save the diagram (choose File-WritePts).  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), sHB ,sSNleft ,sHC , and sSNright.

  1. 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 diagram?
  2. 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?
  3. Repeat for sSNleft < s < sHC, sHC < s < sSNright, and s > sSNright.
D) Superimposing the bifurcation diagram and the numerical solution to the full burster:
  1. 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).
  2. To plot the bifurcation diagram you saved in B), choose GraphicStuff-Freeze-BifDiag.  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.

  1. 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?
  2. 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 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.

  1. 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)?
  2. Repeat with values of 0 < R < 0.5.  Again, what is happening, and can you explain why?
  3. 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?