Most of the figures shown in this chapter were produced by numerically solving the Morris-Lecar equations. We have used a program XPP written by G.B.Ermentrout which uses a variety of numerical integration methods to solve the equations on any UNIX workstation. It is available via anonymous ftp from ftp.math.pitt.edu/pub/bardware and there is an extensive tutorial geared toward neurobiology available on the World Wide Web through GBE's home page http://www.pitt.edu/~phase For the all but the bursting simulations, we have used a fourth order Runge-Kutta algorithm. The bursting simulations employed a variable timestep Gear algorithm [26].
For the Morris-Lecar model, the nullclines can be found explicitly by solving each of the equations for w as a function of v. But this is generally not easy, so that numerical techniques are in order. All of the phase-plane pictures were found by numerically computing the nullclines. This is done by breaking the phase-plane into many small boxes, evaluating the functions on each point and then using a linear interpolation to find the zero contours.
Singular points are found using Newton's method with a numerically computed Jacobian. Once a steady state is found, the Jacobian is computed and the QR algorithm is used to find the eigenvalues. These determine the stability of the singular point.
For certain steady states in the Morris-Lecar model, we want to find special trajectories called the unstable and stable manifolds. This is done by computing an eigenvector for a particular eigenvalue (the eigenvectors are tangent to these manifolds) by inverse iteration. Once the eigenvector is known, the equations are integrated (either forward or backward in time) with initial conditions that are on the eigenvector and slightly off of the singular point.
The qualitative behavior of higher dimensional systems as a parameter is varied can be understood and described compactly by determining the bifurcation diagrams. AUTO [5] was used in this chapter (through its interface with XPP) to trace, essentially automatically, the bifurcation curves as any parameter is varied. This program is able to find all steady states and periodic solutions regardless of their stability. The characteristics of stability, eigenvalues for steady states and Floquet exponents for periodic solutions, are computed along with the frequencies of periodic solutions. Two parameter bifurcation diagrams (similar to Figure 11C) which indicate where steady states and periodic solutions exist and they gain or lose stability can also be computed by AUTO.
The computation of the ``infinitesimal PRC'' is done automatically as part of the XPP package. Essentially, it solves an allied linear equation until a particular periodic orbit is found. The integral required for averaging is automatically computed.
A general and practical reference to many of the above numerical methods, which also leads to more literature is [26].