We have introduced and used some of the basic concepts of the qualititative theory of differential equations to describe the dynamic repertoire of a representative model of excitability. We believe that a geometrical treatment, as in the phase plane, gives one an opportunity to see more clearly and appreciate the underlying qualitative structure of models. One can see which initial conditions, e.g. as the result of a brief perturbing stimulus, will lie in the domain of attraction of any particular stable steady state or limit cycle. This is especially helpful for the design of experiments to switch a multistable system from one mode to another. Analytic methods are also important - for determining and interpreting the stability of solutions (e.g., equation (19) for the Hopf bifurcation) and for approximating aspects of the solution behavior (e.g., equation (25) for phase-locking). Another useful conceptual device is the bifurcation diagram by which we have provided compact descriptions of the system attractors. In several of our illustrations, the bifurcation parameter was I, and in this case the steady-state I-V relation appeared explicitly in the diagram. Any parameter can be used in such considerations such as channel density or synaptic weight as in the Wilson-Cowan network model.
We have shown how a minimal but biophysically reasonable membrane model can be massaged to exhibit robustly a variety of physiologically identifiable firing behaviors. For the simplest two-variable Morris-Lecar model, we illustrated some qualitative differences in threshold behavior when the steady state current-voltage relation is monotonic or N-shaped. In the former case, action potential size may be graded, although generally quite steeply with stimulus strength, and latency for firing is finite. In the latter case, there is a true (saddle point) threshold for action potentials, latency may be arbitrarily long, and intermediate-sized responses are not possible. Correspondingly, for a steady stimulus, the monotonic case leads to onset of oscillations with a well-defined, non-zero frequency (Hopf bifurcation), and with possibly small amplitude (supercritical). In contrast, in the N-shaped case repetitive firing first appears with zero frequency (homoclinic bifurcation). These features are consistent with some of those used by Hodgkin  to distinguish axons with different repetitive firing properties, Class II and Class I, respectively. Additionally, we have provided a geometric interpretation of some common forms of bursting neurons. Many bursters can be dissected into fast dynamics coupled to one or more slow processes that move the fast dynamics between resting and oscillatory states. Coupled and forced oscillators can often be reduced to maps or to continuous low-dimensional systems of phase equations, especially when the interactions are weak.
With regard to implementation of models, there are sophisticated software packages and subroutines available nowadays which help make theoretical experimentation and analysis a real-time endeavor; the ``tooling-up'' time is greatly reduced. Programs like XPP incorporate a mix of numerical integration and analytic formulation (linear stability analysis, bifurcation analysis, averaging - carried out numerically) with graphical representation all in an interactive framework. (The numerical methods employed by XPP are described in Appendix B.) To
set up and get running the Morris-Lecar model should take less than fifteen minutes. Although we have used XPP here for the two-variable model, it can also deal with higher order systems (like the bursting models of Section 4). Stability of steady states can be computed, and time courses plotted interactiviely. The generalization of nullclines to surfaces is not available computationally, but two-variable projections of trajectories from the higher-order phase space can be insightful (e.g., Figure 9). A widely distributed subroutine which we have also found valuable for dissecting nonlinear systems is AUTO  which automatically generates bifurcation diagrams (as in Figures 2,6,7); see Appendix B. The main components of AUTO are included in the package XPP thus making AUTO an interactive program. We used XPP-AUTO for the fast/slow analysis of the bursting models in Section 4. For the evaluation and algebraic manipulation of analytic prescriptions (e.g., lengthy perturbation and bifurcation formulae), many modelers have used symbol manipulation programs like Mathematica and MAPLE with success; see . In the use of numerical packages we advise that one be generally familiar with the methods being employed, and with their limitations. It is not so uncommon to pose a problem which seems to just miss the criteria for suitablity of a given technique - and one should be careful to recognize the symptoms of breakdown of the particular method being used.
Finally, we emphasize the value of using idealized, but biophysically reasonable, models in order to capture the essence of system behavior. If models are more detailed than necessary then identification of critical elements is often obscured by too many possibilities. On the other hand, if justified by adequate biophysical data, more detailed models are valuable for quantitative comparison with experiments. Thus the modeler should be mindful and appreciative of these two different approaches: which one is chosen depends on the types of questions being asked and how much is known about the underlying physiology.