So far, we have looked at a scalar equation for a membrane with a single nonlinear channel. We hinted above that by letting the current vary slowly as a function of the current, we could induce oscillations in the membrane. Real membranes do this but not by directly varing the current. Rather, in addition to the fast inward calcium current, there is an additional voltage controlled outward current due to potassium. That is, the potassium current turns on with high enough voltage and hyperpolarizes the membrane. The equations for this model; the "Morris-Lecar" model are:

where

We will also track the main currents, * ICa * and * IK. *
Click ** HERE ** to see the source.

Let's set up the window so we can explore this system. Time is measured in milliseconds so click (nUmerics) to get the numerics menu. Change (Dt) to .25 and change (Total) to 100 msec. Since the currents can get quite large, change (Bounds) from 100 to 1000. Click (Esc) to return to the main window. Click on (W) (W) to change the dimensions of the window:

- XLo: 0
- XHi: 100
- YLo: -80
- YHi: 120

To understand this, it helps to look at the two currents. Click on
(H) (C) to create another window. Click
(X) and type in ** ICa ** followed by (Enter) and you will see the
large negative calcium current plotted in the window. Now in the
mainwindow, click (G) (A) to add
another curve to the plot. A new window will appear. Type ** IK **
where it asks for the ** Y-axis ** and type in a number from 1-9
for the color (1=red, 9=violet) and then type (Ok) to accept it.
(** Note: ** This is different than
freezing a curve which keeps a ** permanent ** copy of the
curve you have frozen. ** Adding ** a curve just tells XPP to plot
the variables you have added at every simulation. * If you change
parameters or initial conditions, the added curves will change to
reflect this. Frozen curves do not change.*)

You will not see the new curve since XPP has not given you the best window. You can manually increase the window, or let XPP do it. Click (W) (F) to fit the window. (NOTE: This only works in versions of XPP higher than 1.6.) Otherwise just use (W) (W) and make the vertical scales from -380 to 270.

Click (T) for text and at the prompt type in "ICa" and (Enter). Using the mouse, place the text next to the big negative calcium current and click it. Repeat this for "IK." If you have done this, you will see something like:

Now we can see what * excitability * means. Once the voltage
crosses threshold there is a large inward calcium current. After this
reaches its peak causing the voltage to increase, the large outward
potassium current comes into play. This is sufficient to cause the
membrane to repolarize back to rest.

Click inside the big graphics window showing the voltage. Click on (H)
(K) to kill all extra windows and answer (Yes). This will get rid of
the little window. Now click on (U) to go into the (nUmerics) window
and change the (Total) time to 250 msec and hit (Esc) to return to the
main window. Click (V) (2) to get a Two-d
view of the voltage and the potassium gating variable, * w
*. A new window will appear. You must tell XPP the variables for
the X and Y axes, the (optional) labels and dimensions of the window
(although this can be done with the (Window) command.) We want the
"X-axis" to have voltage and the "Y-axis" to have the gating
variable. Fill it in as:

- X-axis: V
- Y-axis: w
- XMin: -80
- YMin: -.25
- XMax: 120
- YMax: 1.25
- XLabel: V
- YLabel: w

This is called a ** phase-plane ** picture and it is what the "PP"
in XPP stands for. The ** phase-space ** for the Morris-Lecar
model is two-dimensional and this provides a nice compact way of
seeing what the solutions do as a function of time.
Click (I) (N) and set V(0)=-20 and leave w as is. Note the small
trajectory thet goes to the left. Now initialize V(0)=-10 and see the
difference. The trajectory loops through the two-dimensional plane
before returning to rest. Note that the phase-plane view shows you
that the voltage increases and the gating variable slowly increases.
The voltage reaches a maximum but the gating variable continues to
increase reaching a maximum and then decreasing. The voltage
udershoots the equilibrium before returning to rest.

Now click on (S) (M) to use the mouse to make an initial guess for a singular point or equilibrium. Click the mouse near where the two trajectories have ended. Answer (No) to the "Print Eigenvalues" question. A small circle will appear in the phase-plane that shows the position of the equilibrium. The stability window tells you you have a rest state at V=-60.89 and w=.014, that it is stable, and that the two eigenvalues are complex.

Click on (I) (M) to choose initial conditions in the phase-plane with the mouse. Click anywhere in the plane and the trajectory will be plotted. In particular, try a bunch of initial conditions near w=0 along the V-axis. Turn off the Bell by clicking (F) (B). Now click (I) (L) to get a trajectory that starts close to the initial condition. ((I) (L) means use the (Last) value of the previous trajectory as the starting point of this one.) Now click (I) (R) to do a range integration. A new window will appear. Choose V to be the variable to range over. Choose 10 Steps, with the start at -20 and the end -10. Answer yes to all three questions and click (Ok). You will see 10 colored trajectories with a clear point indicating threshold. You should see something like this:

(* If this is nothing like what you see, try going through the
instructions again starting at the beginning. *)

To see what is causing the threshold behavior for the neuron, it is
very helpful to look at the **
nullclines ** which are curves where
* dV/dt=0 * and * dw/dt=0.* The intersection of two
nullclins is an equilibrium point. Click (N)(N) to make new
nullclines. The "V" nullcline where V does not change is RED and the
"w" nullcline, where w does not change is green. They intersect in one
point, the equilibrium. Their meaning can be discerned from the
figure below:

The nullclines are labeled and divide the phase-plane into 4 regions
according to whether the voltage V and gating variable, w,
are increasing or decreasing. THis lets you qualitatively see what
happens. Use the (I) (M) command to start with initial conditions in
each of the 4 regions of the phase-plane. Click on (E) to erase the
screen and then on (N) (R)
to restore the nullclines without having to
recompute them. Now use the mouse to choose initial conditions along
the line * w=0 * (the x-axis) near the point where the
V-nullcline intersects the x-axis. Try some intial data to the left
and to the right of the curve. * The nullcline acts almost like a
threshold. *

Now set I=120 by typing (P) and then I and then the value for it. Hit
(Enter) twice. If you have the
parameter window open then just click on I in the window and tyoe
the new value in the command line.
Once again use the mouse to set the initial data. Click (I) (M) and
then click the mouse in the plane. You should see a
closed circle in the phase plane, a ** limit cycle ** for which
both variables are periodic. Open another window and plot V as a
function of t. Click (Halfwindow) (Create). Then click (Xivst) and
type in V followed by (Enter). You should se a nice oscillation of
the membrane potential. Click on (Xivst) and type in * w * and
see that it oscillates as well.

- Look at ICa and IK in a two-dimensional phase-space.
- Look at the
nullclines in the
**V-w**plane when**I=120**and note that they intersect on the part of the**V-nullcline**with a positive slope. Find the stability of the singular point where the nullclines intersect. What is the nature of the eigenvalues? Change**I=80**and look at the nullclines and singular point again. How has the rest state changed? What behavior occurs for values of the parameter at which the rest state is*unstable?* - Try to find the value of
**I**at which the singular point goes from stable to unstable. (*Hint:*Once you have found an equilibrium point using the (S) (M) option, try using (S) (G) after you change**I**slightly. )

We will use AUTO to find that branch of periodics as well as discover
some other interesting behavior. Set ** I=0 ** and initialize **
V=-60.899, w=.01487.** Now click (F) (A) to invoke AUTO. Now, setr
up the axes by clicking on (Axes), choosing (Hi-Lo) and filling in the
window as:

- Y-axis: V
- Main Parm: i
- 2nd Parm: vl
- Xmin: -10
- Ymin: -80
- Xmax: 250
- Ymax: 120

- Ds: 0.02
- Dsmin: 0.0001
- Dsmax: 5
- Par Min: -10
- Par Max: 245
- Norm Max: 1000

We are almost ready to roll!!. AUTO and XPP let us record aspects of
the bifurcation curves at chosen values of parameters and periods of
limit cycles. We want to get information at ** I=100 ** and,
say, ** I=120.**
This is done by clicking on (Usr Period). Do this and
choose 2. A window pops up. Fill in the first two entries like:

- Uzr1: i=100
- Uzr2: i=120

Now we can run! Click on (Run) and choose (Steady state). A thick line which is thin in the middle is drawn across the screen. The thick parts are where the steady state is stable, the thin, unstable.

How does the steady state lose stability? Click on (Grab) and
move the cross-hairs until they are close to the instability. Look at
the information line as you do this. You should see a symbol ** HB
** at around ** I=101.8 ** How does this compare with your
estimate of the value of current which leads to instability
in the homework? The symbol ** HB ** means that there is a ** Hopf
bifurcation ** so we can expect to see a branch of solutions
emanating from the rest state. When the cross-hair is on the ** HB
** point, click
(Enter) to accept it. When you accept a point in AUTO, this sets the
parameters and dependent variables to the initial condition
corresponding to the point. It also gives AUTO some necessary
information.

Now that we have grabbed the Hopf point, we can try to find the
periodic orbit. * Note: The next command may induce AUTO to
loop around so once it looks like AUTO is retracing stuff it already
found, click on (ABORT). * So forewarned, click on (Run). Choose
(Periodic) to get AUTO to trace the orbit. Stop it if it looks like
its retracing! When the diagram is complete it will look like:

Filled circles represent stable periodic orbits and open circles are
unstable. Notice how the unstable branch that emanates from the Hopf
bifurcation bends to the left and thet the steady state solution loses
stability to the right. If the direction of the bifurcating branch is
opposite the direction at which stability of the main branch is lost,
the branch is called ** subcritical ** and it is generally **
unstable. ** Note that there are two subcritical Hopf bifurcations;
the second one is for a high current of about 220. (It is subcritical
because the steady state branch loses the stability going * left
* and the periodic branch goes * right. *)

Click (Grab). Using the (Tab) key, slowly click through the special points. Note the labels `
UZ ` and ` LP `. The label ` UZ ` is a user-defined
point and since we picked **I=100 ** and ** I=120 **, it means
that ** I ** must be one of these two values. The label ` LP
` means that a limit-point or saddle node has been found. Since
it occurs on a branch of periodic solutions, that means that two limit
cycle orbits coalesce at this point. For this
particular model, the left most limit point occurs at ** I =
95.72. **

In many biological examples, subcritical branches from Hopf
bifurcations * turn around. * The Morris-Lecar model, here, is
such an example.
Focus your attention on the region between the Hopf bifurcation and
the left-most limit point, ** 95.72 < I < 101.8 ** For
currents in this interval, there are 3 different possible
solutions. (Imagine a vertical line drawn in this interval; it passes
through three branches.) There is a * stable steady state, unstable
periodic, stable periodic. *

- Using (Grab) load the second point on
the bifurcation diagram for which
**I=100**(label 8). - Now click (File) (Import Orbit). Special points have more information than regular points. Grabbing any point changes the parameter and initial conditions in XPP. But you can actually load the complete trajectory of a periodic orbit if it is a special one.
- In the XPP main window, you should have a phase-plane projection
of V and w. If not, then do so now
with
*-80 < V < 120*and*-.25 < w < 1.25.*Click on (Restore) in the XPP main window. (Since the nullclines are probably wrong, erase the screen by clicking on (Erase) then draw new nullclines for this value of I by clicking (N) (N), then click (R) to restore.) You should see a very small periodic orbit surrounding the fixed point marked by the intersection of the nullclines. Check the stability of the rest state by clicking (S) (M) in the XPP window and clicking the mouse at the intersection -- answer (No) to the eigenvalue question. It is stable. - Lets get a permanent copy of the small orbit. Click (G) (F) (F) to freeze the small orbit. Choose color 4 and make the key have the label "UPO."
- Go back to the AUTO window and click grab again. This time load the third instance for which the current is 100. This is labeled point 10. Again, click (F) (I) in the AUTO window to import this orbit.
- Go back to the main XPP window and freeze this using (G) (F) (F) and choosing color 7 and labeling the key, "SPO". You should see two periodics in the picture, one stable and one unstable.
- Finally, choose a variety of initial conditions by clicking (I) (M) in the XPP window and clicking the mouse inside the small periodic, between the two periodics, and outside the large periodic. (It might be easier to see if you zoom in so that the picture contains the two periodics. Click (W) (Z) and use the mouse key to zoom in. )

If you have followed these instructions, then you will have
discoverd another type of ** bistability ** between a stable fixed
point and a stable periodic orbit. All initial conditions * inside
* the unstable periodic orbit end up at the stable fixed point.
All initial data * outside * the unstable periodic end up on the stable
periodic state. The unstable periodic orbit is a threshold for
repetitive firing. Create another window and plot V as a function of
t for some represetative trajectories.

This concludes the tutorial on periodic, bistability, and AUTO. You can play around some more with this diagram (for example, pick some currents at the high end where there is also bistability.) When you are done, reset the bifurcation diagram to delete unwanted files from the disk: In the AUTO window, reset the diagram by clicking (File) (Reset) and answering (Yes). Exit XPP or try the next homework which is quite extensive.

- V3=12
- V4=17.4
- I=0
- phi=.066666667

- Dt=.25
- Bound=1000
- Total=250

- Look at the nullclines for I=0, I=30, I=50. Look at trajectories and determine the stability of all fixed points. At what value of I does repetitive firing occur? For I=30, what determines the threshold? Try the following with I=30. Draw the nullclines. Click on (S) (M) to use the mouse to find singular points and choose the middle intersection of the nullclines. Answer (No) to the eigenvalue question and (Yes) to the invariant sets. Watch the trajectories that appear. (Click Esc to make it draw the next one, -- it will try to draw 4 of them) Note the pair that come out almost vertically. Click some initial conditions on either side of them and note the behavior. They determine the threshold.
- Now, set I=0 and initialize V to -59.474 and w to 0. Fire up the AUTO window if it is not already there (F) (A). If it is there,make sure the diagram has been reset (Click (File) (Reset Diagram)). Now set up the axes and numerics just like the example above. Run AUTO starting from a steady state. How is the steady state I-V curve different from the first example? Notice that even though it is "kinked" there are never two stable rest states for the same parameters.
- How does the upper steady state gain/lose stability? Use (Grab) to find out if there is a Hopf bifurcation. Grab it and trace out the curve of periodics. Compare this scenario with the above one. Look at the frequency as a function of I. Concentrate on the stable branch. What is the minimal frequency?
- In this scenario, the periodic orbit emerges with large amplitude and zero frequency. It is called Type I behavior. The limit cycle disappears exactly at the value of current for which the two lower steady states come together.

- Table of Contents
- Main Menu Items
- ODE Files and Examples
- Numerics Menu
- File Menu
- Freeze Menu
- AUTO Menu
- Data Browser
- I/O and Hardcopy
- XPP Basics
- Nonlinear ODEs
- Two-dimensions
- Three-dimensions and Beyond
- Phase Equations
- Chaos