%\renewcommand{\here}{/a2.encfigs/}
\renewcommand{\here}{}
\setcounter{chapter}{5}

\chapter{Model Neurons}

\section{Levels of Neuron Modeling}

A great deal is known about the biophysical mechanisms responsible for
generating neuronal activity, and the mathematical descriptions of
these mechanisms presented in chapter 5 provide a basis for
constructing neuron models. These models range from highly detailed
descriptions involving thousands of coupled differential equations to
greatly simplified caricatures useful for large network modeling.  In
this chapter, we present a number of different approaches to modeling
neuronal activity that vary considerably in their level of detail.
Choosing the most appropriate modeling level for a given research
problem requires a careful assessment of the experimental information
available and a clear understanding of the research goals.
Oversimplified models can, of course, give misleading results, but
excessively detailed models can obscure interesting results beneath
inessential and unconstrained complexity.

In modeling neurons, we must deal with two basic types of complexity:
elaborate morphologies and conductances with rich dynamics.  As
discussed in chapter 5, the membrane potential of a neuron with a
complex morphology depends on both time and spatial position, $V(x,
t)$. The cable equation that determines $V(x, t)$ was discussed in
chapter 5, but the analytic methods presented there can only be
applied when the membrane conductances are voltage-independent.  When
the complexities of real membrane conductances are included, the
membrane potential must be computed numerically.  This is done by
breaking up the neuron into separate regions or compartments and
approximating the continuous membrane potential $V(x, t)$ by a
discrete set of values representing the potentials within the
different compartments. This assumes that each compartment is small
enough so that there is negligible variation in the membrane potential
across it. The precision of such a multi-compartmental description
depends on the number of compartments used and on the compactness of
the neuron being modeled.  Figure \ref{B2:compneuron} shows a
schematic diagram of a cortical pyramidal neuron, along with a series
of compartmental approximations of its structure.  The number of
compartments used can range from thousands in some models to one for
the description at the extreme right of figure \ref{B2:compneuron}.
In single compartment models, the entire neuron is described by a
single membrane potential and all position-dependent variations are
ignored.  This is only valid if the neuron is fairly compact, or if
the behavior being studied is not strongly affected by potential
differences between different parts of the neuron.

\begin{figure}[hhtt!]
\psfig{figure=figs/compneuron.eps,width=4.5in,clip} 
\caption{A sequence of approximations for the structure
    of a neuron.  The neuron is represented by a variable number of
    discrete compartments each representing a region that is described
    by a single membrane potential. The connectors between
    compartments represent resistive couplings.  The simplest
    description is the single-compartment model furthest to the right.
    (Neuron diagram from Haberly, 1990)}
\label{B2:compneuron}
\end{figure}

One of the most problematic features in neuron modeling is the
tremendous dynamic range of neuronal membrane and synaptic
conductances.  The fastest conductances in a typical neuron, those
responsible for initiating the action potential, activate in less than
a millisecond.  Slower conductances have time constants in the range
of hundreds of milliseconds.  Various modulatory processes can affect
conductances over times varying from minutes to hours, or even longer.
Thus, the complete dynamic range extends over several orders of
magnitude.  In most models, this is collapsed to a much smaller range
to keep the time needed to numerically integrate the equations of the
model reasonable.

The ionic conductances that give a neuron its electrical properties
are time-dependent nonlinear functions of voltage and other factors.
Models that treat these aspects of ionic conductances are known as
conductance-based models.  The basic formalism developed by Hodgkin
and Huxley to describe the \Na\ and \K\ conductances responsible for
generating action potentials (discussed in chapter 5) is used to
represent most of the additional conductances encountered in neuron
modeling.  Neuron models based on these mathematical descriptions of
membrane conductances can reproduce the rich and complex dynamics of
real neurons quite accurately.  In this chapter, we discuss both
single- and multi-compartment conductance-based models.

When large networks of neurons are modeled, simplified descriptions of
the individual neurons are often employed.  This speeds up simulations
of large networks, simplifies the analysis of their activity patterns,
and avoids the problem of setting large numbers of free parameters
that are not constrained by experimental measurements.  We discuss a
particularly useful simplified model of a spiking neuron, the
integrate-and-fire model, quite extensively.  We conclude the chapter
by analyzing a complementary approach, neuron models based on
firing-rate descriptions that provide the basis for the network models
discussed in chapters 7 and 8\@.

\section{Single-Compartment Models}

In any neuron model, the rate of change of the membrane potential
within a portion of the neuron is determined by the total current
flowing into that region (see chapter 5).  For a single-compartment
model, the region in question is the entire neuron, and the relevant
currents are those arising from all the membrane and synaptic
conductances of the cell.  The sum of all these currents is set equal
to the total capacitance of the neuron times the rate of change of the
membrane potential. Because membrane and synaptic currents are usually
expressed as currents per unit area of membrane, it is more convenient
to re-express this relation so that the total current per unit area of
membrane is equal to the capacitance per unit area, $c_{\m}$ (usually
taken to be 10 nF/mm$^2$), times the rate of change of the membrane
potential.  We also include a current $I_{\emy}$ coming from an
external electrode inserted into the neuron.  Electrode currents are
not typically expressed as current per unit area of membrane, so to
match the other currents we must divide the injected current by the
total surface area of the neuron, $A$\@.  Putting all this together,
the basic equation for all single-compartment models is
\begin{equation}
c_{\m}\frac{dV}{dt} = -i_{\m} - i_{\s}+ \frac{I_{\emy}}{A}
\label{B2:onec}
\end{equation}
where $i_{\m}$ is the total current per unit area from all membrane
conductances, and $i_{\s}$ is the total current per unit area due to
synaptic conductances.  By convention, current that enters the neuron
through an electrode is defined as positive-inward, whereas membrane
and synaptic currents are defined as positive-outward.  This explains
the different signs for the currents in equation \ref{B2:onec}.
Single-compartment neuron models all use equation \ref{B2:onec}, but
they differ in the degree of detail with which the dynamics of the
membrane and synaptic conductances are treated.

\subsection*{Conductance-Based Models} 

The formalism developed by Hodgkin and Huxley (1952) (see chapter 5)
to describe the fast \Na\ and delayed-rectifier \K\ conductances in
the squid giant axon has, by now, been applied successfully to dozens
of conductances in a wide variety of preparations.  On the basis of
these conductance descriptions, it is possible to construct complex
and reasonably accurate models of neurons displaying a wide range of
activity patterns.  Efficient numerical methods make it possible to
run computer simulations that reproduce the activity of real neurons
rapidly and with a high degree of accuracy.

\begin{figure}[hhtt!]
\epsfig{figure=figs/1comp.eps, width=4.5in, clip}
\caption{\small{The equivalent circuit for a one-compartment
    conductance-based neuron model.  The neuron is modeled, at left,
    as a single compartment with a synapse and a current injecting
    electrode. At right is the equivalent circuit. The ellipses stand
    for an unspecified number of membrane conductances. The circled s
    indicates a synaptic conductance that depends on the activity of a
    presynaptic neuron.  A single synaptic conductance is indicated,
    but, in general, there may be several different types. The circled
    V indicates a voltage-dependent conductance.  $I_{\emy}$ is the
    current passing through the electrode, and $A$ is the surface area
    of the neuron.}} \label{B2:1comp}
\end{figure}
The membrane potential for a single-compartment conductance-based
neuron model is determined by integrating equation \ref{B2:onec} (see
Appendix A).  The membrane and synaptic currents per unit area are
written in terms of the membrane and synaptic conductances, open
probabilities, or gating variables (see chapter 5) as
\begin{equation}
i_{\m} + i_{\s} = \sum_i g_i(V-E_i) = \sum_i\overline g_iP_i(V - E_i)  
= \sum_i\overline g_im_i^{p_i}h_i^{q_i}(V - E_i)
\label{B2:memcurr}
\end{equation}
where the parameters $\overline g_i$ are fixed maximal conductances
(the conductance per unit area when $P=1$).  Note that we have
included both membrane and synaptic conductances in this sum, because
both are describe similarly in terms of gating variables and reversal
potentials.  One of the membrane conductances is the leakage
conductance with $P=1$.  Figure \ref{B2:1comp} shows the equivalent
circuit for a generic one-compartment model with a number of different
conductances.

The dynamic variables of a single-compartment conductance-based model
are the membrane potential, $V$, and the activation and inactivation
gating variables for the different conductances.  We have
distinguished the different gating variables in equation
\ref{B2:memcurr} by the subscript $i$, but often they are
distinguished instead by using letters other than $m$ and $h$.  For
example, the activation variable for the delayed-rectifier \K\ 
conductance is conventionally denoted by $n$ (see chapter 5).  As
discussed in chapter 5, all of the gating variables (whether denoted
by $m$, $h$, or some other letter) are determined by equations of the
form
\begin{equation}
\tau_z(V)\frac{dz}{dt} = z_\infty(V) - z
\label{B2:zeqn}
\end{equation} 
where we have used the letter $z$ to denote a generic gating variable.
The functions $\tau_z(V)$ and $z_\infty(V)$ are determined from
experimental data. For some conductances, these are written in terms
of the open and closing rates $\alpha_z(V)$ and $\beta_z(V)$ (see
chapter 5) as
\begin{equation}
  \tau_z(V) =\frac{1}{\alpha_z(V)+\beta_z(V)}
  \;\;\;\;\;\mbox{and}\;\;\;\;\; z_\infty(V) =
  \frac{\alpha_z(V)}{\alpha_z(V)+\beta_z(V)} \, .
\label{B1:alphabeta}
\end{equation}
We have written $\tau_z(V)$ and $z_\infty(V)$ as functions of the
membrane potential, but for \Ca -dependent currents they also depend
on the internal \Ca\ concentration.  For synaptic conductances, the
gating variables depend on the concentration of transmitter in the
synaptic cleft, but they are often written directly in terms of the
presynaptic membrane potential or the times when presynaptic action
potentials occur (see chapter 5).  A method for numerically
integrating the equations for the gating variables of a
conductance-based model is described in Appendix B\@.

In the following sections, some basic features of conductance-based
models are presented in a sequence of examples of increasing
complexity.  We do this to illustrate the effects of various
conductances and combinations of conductances on neuronal activity.

\subsubsection*{The Connor-Stevens Model}

The Hodgkin-Huxley model of action potential generation, discussed in
chapter 5, was developed on the basis of data from the giant axon of
the squid, and we present a multi-compartment simulation of action
potential propagation using this model in a later section.  The
Connor-Stevens model provides an alternative description of action
potential generation.  Like the Hodgkin-Huxley model, it contains fast
\Na , delayed-rectifier \K , and leakage conductances.  The fast \Na
and delayed-rectifier \K\ conductances have somewhat different
properties from those of the Hodgkin-Huxley model, in particular
faster kinetics, so the action potentials are briefer.  In addition,
the Connor-Stevens model contains an extra \K\ conductance, called the
A-current, that is transient.

The membrane current in the Connor-Stevens model is
\begin{equation} 
  i_{\m} ={\overline g}_L(V-E_{\Lmy})+{\overline
    g}_{\Namy}m^3h(V-E_{\Namy}) +{\overline g}_Kn^4(V-E_{\Kmy})+
  {\overline g}_{\A} a^3b(V-E_{\A})
\end{equation} 
where ${\overline g}_{\Lmy} = 3$ $\mu$S/mm$^2$ and $E_{\Lmy}$ = -17 mV
are the maximal conductance and reversal potential for the leak
conductance, and ${\overline g}_{\Kmy} = 200$ $\mu$S/mm$^2$,
${\overline g}_{\Namy} = 1200$ $\mu$S/mm$^2$, ${\overline g}_{\A} =
477$ $\mu$S/mm$^2$, $E_{\Kmy}$ = -72 mV, $E_{\A}$ = -75 mV and
$E_{\Namy}$ = 55 mV are similar parameters for the active
conductances\@.  The gating variables, $m$, $h$, $n$, $a$, and $b$,
are determined by equations of the form \ref{B2:zeqn}.  The $\alpha$
and $\beta$ functions for the \Na\ and delayed-rectifier \K\ 
conductances are (with $\alpha$ and $\beta$ in units of 1/ms and $V$
in units of mV)
\begin{eqnarray}
  \!\!\!\!\!\!\!\!\!\!\alpha_m=\frac{0.38(V+29.7)}{1-\exp[-0.1(V+29.7)]}&
  \beta_m=15.2\exp[-0.0556(V+54.7)] \nonumber
\\
\!\!\!\!\!\!\!\!\!\!\alpha_h=0.266\exp[-0.05(V+48)]&\beta_h= 3.8/(1+\exp[-0.1(V+18)])
\nonumber
\\ 
\!\!\!\!\!\!\!\!\!\!\alpha_n=\frac{0.02(V+45.7)}{1-\exp[-.1(V+45.7)]}&
\beta_n=0.25\exp[-0.0125(V+55.7)] \, .
\end{eqnarray}
The A-current is described directly in terms of the asymptotic values
and $\tau$ functions for its gating variables (with $\tau_a$ and
$\tau_b$ in units of ms and $V$ in units of mV),
\begin{equation}
a_\infty=\left[\frac{0.0761\exp[.0314(V+94.22)]}{1+\exp[0.0346(V+1.17)]}\right]^{1/3}
\end{equation}
\begin{equation} 
\tau_a=0.3632 + 1.158/(1+\exp[0.0497(V+55.96)])
\end{equation}
\begin{equation}
b_\infty=\left[\frac{1}{1+\exp[0.0688(V+53.3)]}\right]^4
\end{equation}
and
\begin{equation} 
\tau_b=1.24 + 2.678/(1+\exp[0.0624(V+50)]) \, .
\end{equation}

\begin{figure}
\epsfig{figure=figs/csRates.eps, width=4.5in}
\caption{\small{Firing of action potentials in the Connor-Stevens
    model.  A) Firing rate as 
    a function of injected current.  The firing rate rises
    continuously from zero as the current increases beyond the
    threshold value.  B) An example of action potentials generated by
    constant current injection.  C) Firing rate as a function of
    injected current when the A-current is turned off.  The firing
    rate now rises discontinuously from zero as the current increases
    beyond the threshold value.  D) Delayed firing due to
    hyperpolarization.  The neuron was held hyperpolarized for a
    prolonged period by injection of negative current.  At $t = 50$
    ms, the negative injected current was switched to a positive
    value.  The A-current delays the occurance of the first action
    potential.}}
\label{B2:csRates}
\end{figure} 
Figure \ref{B2:csRates} illustrates action potential generation in the
Connor-Stevens model.  In the absence of injected external current or
synaptic input, the membrane potential of the model remains constant
at a resting value of $-68$ mV\@.  With constant positive current
injection ($I_{\emy}> 0$), the model neuron depolarizes until a
threshold current is reached and the model generates action
potentials.  Figure \ref{B2:csRates}A shows how the firing rate of the
model depends on the magnitude of the injected current relative to the
threshold value.  The firing rate rises continuously from zero and
then increases roughly linearly for currents over the range shown.
For larger currents, the firing rate increases more slowly and
ultimately saturates.  Figure \ref{B2:csRates}B shows an example of
action potential generation due to constant current injection.

The fast \Na\ and delayed-rectifier \K\ conductances in the
Connor-Stevens model generate action potentials in the same way they
do in the Hodgkin-Huxley model (see chapter 5).  What is the role of
the addition A-current?  Figure \ref{B2:csRates}C shows the firing
rate as a function of injected current for the Connor-Stevens model
with the maximal conductance of the A-current set to zero.  The
leakage conductance and reversal potential have been adjusted to keep
the same resting potential and input resistance as in the original
model.  The firing rate is clearly much faster with the A-current
turned off.  In addition, the transition from no firing for currents
less than the threshold value to firing with suprathreshold currents
is different when the A-current is eliminated.  Without the A-current,
the firing rate jumps discontinuously to a nonzero value rather than
rising continuously.  Neurons with firing rates that rise continuously
from zero as a function of injected current are called type I, and
those with discontinuous jumps in their firing rates at threshold are
called type II\@.  An A-current is not required to produce a type I
response but, as figures \ref{B2:csRates}A and \ref{B2:csRates}C show,
it plays this role in the Connor-Stevens model. The Hodgkin-Huxley
model produces a type II response.


Another effect of the A-current is illustrated in figure
\ref{B2:csRates}D\@.  Here the model neuron was held hyperpolarized by
negative current injection for an extended period of time, and then
the current was switched to a positive value.  While the neuron is
hyperpolarized, the A-current deinactivates; that is, the variable $b$
increases toward one.  When the injected current switches sign and the
neuron depolarizes, the A-current first activates and then
inactivates.  This delayes the first spike following the change in the
injected current.

\subsubsection*{Postinhibitory Rebound and Bursting}

The range of responses exhibited by the Connor-Stevens model neuron
can be extended by including a transient \Ca\ conductance.  The
conductance we use was modeled by Huguenard and McCormick (1992) on
the basis of data from thalamic relay cells.  The membrane current due
to the transient \Ca\ conductance is expressed as
\begin{equation}
i_{\CaT} = \overline g_{\CaT}M^2H(V - E_{\Camy})
\end{equation}
with, for the example given here, $\overline g_{\CaT} = 13$
$\mu$S/mm$^2$ and $E_{\Camy} = 120$ mV\@. The gating variables for the
transient \Ca\ conductance are determined by (with $\tau_M$ and
$\tau_H$ in ms and $V$ in mV)
\begin{equation}
M_{\infty} = \frac{1}{1 + \exp\left(-(V + 57)/6.2\right)}
\end{equation}
\begin{equation}
H_{\infty} = \frac{1}{1 + \exp\left((V + 81)/4\right)}
\end{equation}
\begin{equation}
\tau_M = 0.612 + \left(\exp\left(-(V + 132)/16.7\right) + 
\exp\left((V + 16.8)/18.2)\right)\right)^{-1}
\end{equation}
and
\begin{equation}
\tau_H = \left\{\begin{array}{ll}
\exp\left((V + 467)/66.6\right) & \mbox{if $V < -80$ mV}\\
28 + \exp\left(-(V + 22)/10.5\right)  & \mbox{if $V \geq -80$ mV} \, .
\end{array} \right.
\end{equation}
\begin{figure}
\epsfig{figure=figs/caSpike.eps, width=3.5in}
\caption{\small{A burst of action potentials due to rebound from hyperpolarization.  The
    model neuron was held hyperpolarized for an extended period (until
    the conductances came to equilibrium) by injection of constant
    negative external current.  At $t= 50$ ms, the injected current
    was set to zero and a burst of \Na\ spikes was generated due to an
    underlying \Ca\ spike.}}
\label{B2:caSpike}
\end{figure} 

A transient \Ca\ conductances acts, in many ways, like a slower
version of the transient \Na\ conductance that generates action
potentials.  Instead of producing an action potential, the transient
\Ca\ conductance generates a slower transient depolarization sometimes
called a \Ca\ spike.  This transient depolarization causes the neuron
to fire a burst of action potentials, which are \Na\ spikes riding on
the slower \Ca\ spike.  Figure \ref{B2:caSpike} shows such a burst and
illustrates one way to produce it.  In this example, the model neuron
was hyperpolarized for an extended period and then released from
hyperpolarization by setting the current to zero.  During the
prolonged hyperpolarization, the transient \Ca\ conductance
deinactivated.  When the injected current was set to zero, the
resulting depolarization activated the transient \Ca\ conductance and
generated a burst of action potentials.  The burst in figure
\ref{B2:caSpike} is delayed due to the presence of the A-current in
the original Connor-Stevens model, and it terminates when the \Ca\ 
conductances inactivates.  Generation of action potentials in response
to release from hyperpolarization is called postinhibitory rebound
because, in a natural setting, the hyperpolarization would be caused
by inhibitory synaptic input, not by current injection.

The transient \Ca\ current is an important component of models of
thalamic relay neurons.  These neurons exhibit different firing
patterns in sleep and wakeful states.  Action potentials tend to
appear in bursts during sleep. Figure \ref{B2:xjmod} shows an example
of three states of activity of a model thalamic relay cell due to Wang
(1994) that has, in addition to fast \Na , delayed-rectifier \K , and
transient \Ca\ conductances, a hyperpolarization activated
mixed-cation conductance and a persistent \Na\ conductance. The model
is silent or fires action potentials in a regular pattern or in bursts
depending on the level of current injection.  In particular, injection
of small amounts of negative current leads to bursting.  This occurs
because the hyperpolarization due to the current injection
deinactivates the transient \Ca\ current and activates the
hyperpolarization activated current.
\begin{figure}
\epsfig{figure=figs/xjmod.eps, width=4.5in}
\caption{\small{Three activity modes of a model thalamic neuron.  Upper panel: with no
    injected current the model is silent.  Middle panel: when a
    positive current is injected into the model neuron, it fires
    action potentials in a regular periodic pattern. Lower panel: when
    negative current is injected into the model neuron, it fires
    action potentials in periodic bursts.  (Adapted from Wang, 1994)}}
\label{B2:xjmod}
\end{figure} 

Neurons can fire action potentials either at a steady rate or in
bursts even in the absence of current injection or synaptic input.
Periodic bursting is a common feature of neurons in central patterns
generators, which are neural circuits that produce periodic patterns
of activity to drive rhythmic motor behaviors such as walking,
running, or chewing.  To illustrate periodic bursting, we consider a
model constructed to match the activity of neurons in the crustacean
stomatogastric ganglion (STG), a neuronal circuit that controls
chewing and digestive rhythms in the foregut of lobsters and crabs
(this is a variant of the model of Turrigiano, LeMasson, and Marder,
1995 due to Z. Liu and M. Goldman).  The model contains fast \Na ,
delayed-rectifier \K , A-type \K , and transient \Ca\ conductances
similar to those discussed above, although the formulae and parameters
used are somewhat different.  In addition, the model has a
\Ca-dependent \K\ conductance.  Due to the complexity of the model, we
do not provide complete descriptions of its conductances except for
the \Ca-dependent \K\ conductance, a type of conductance we have not
yet discussed.  This conductance plays a particularly significant role
in the model.

The \Ca -dependent \K\ current is given by
\begin{equation}
i_{\KCa} = \overline g_{\KCa}c^4(V - E_{\Kmy})
\end{equation}
where (with $V$ in mV)
\begin{equation}
c_\infty = \left(\frac{[\mbox{\Ca}]}{[\mbox{\Ca}] + 3 \mu\mbox{M}}\right)
\frac{1}{1 + \exp(-(V + 28.3)/12.6)}
\end{equation}
and (with $\tau_c$ in ms and $V$ in mV)
\begin{equation}
\tau_c = 90.3 - \frac{75.1}{1 + \exp(-(V + 46)/22.7)} \, .
\end{equation}
Note that $c_\infty$ depends on both the membrane potential and the intracellular \Ca\
concentration, $[\mbox{\Ca}]$. 

The intracellular \Ca\ concentration is computed in this model using a
greatly simplified description in which rises in intracellular \Ca\ 
are caused by influx through membrane \Ca\ channels, and \Ca\ removal
is described by an exponential process.  The resulting equation for
the intracellular \Ca\ concentration, [\Ca ], is
\begin{equation}
\frac{d[\mbox{\Ca}]}{dt} = -\gamma i_{\Camy} - \frac{[\mbox{\Ca}]}{\tau_{\Camy}} \, .
\end{equation}
Here $i_{\Camy}$ is the total \Ca\ current per unit area of membrane, $\tau_{\Camy}$ is the
time constant determining the \Ca\ uptake rate, and 
\begin{equation}
\gamma = \left(\frac{A}{\mbox{Vol}F}\right)\left(\frac{10^6
\mbox{ mm}^3}{\mbox{liter}}\right) 
\end{equation}
where $A$ is the total surface area of the cell, Vol is its total
volume, and $F$ is the Faraday constant.  The factor $\gamma$ converts
the \Ca\ current, which is measured as a current per unit area, into a
volume flux of \Ca\ measured in mol units.

\begin{figure}
\epsfig{figure=figs/burster.eps, width=4.5in}
\caption{\small{Periodic bursting in a model of a crustacean stomatogastric ganglion neuron. 
    From the top, the panels show the membrane potential, the \Ca\ 
    conductance, the intracellular \Ca\ concentration, and the \Ca
    -dependent \K\ conductance.  The \Ca -dependent \K\ conductance is
    shown at an expanded scale so the reduction of the conductance due
    to the falling intracellular \Ca\ concentration during the
    interburst intervals can be seen.  In this example, $\tau_{\Camy}
    = 200$ ms. (Simulation by M. Goldman.)}}
\label{B2:burster}
\end{figure} 
Figure \ref{B2:burster} shows the model firing action potentials in
bursts.  As in the models of figures \ref{B2:caSpike} and
\ref{B2:xjmod}, the bursts are transient \Ca\ spikes with action
potentials riding on top of them.  The \Ca\ current during these
bursts causes a dramatic increase in the intracellular \Ca\ 
concentration.  This activates the \Ca -dependent \K\ current which,
along with the inactivation of the \Ca\ current, terminates the burst.
The interburst interval is determined primarily by the time it takes
for the intracellular \Ca\ concentration to return to a low value,
which deactivates the \Ca -dependent \K\ current and allows another
burst to be generated.  Although figure \ref{B2:burster} shows that
the conductance of the \Ca -dependent \K\ current reaches a low value
immediately after each burst, this initial dip is too early for
another burst to be generated at that point in the cycle.

\subsubsection*{Activity-Dependent Conductances}

Conductance-based neuron modeling is an area of computational
neuroscience where experimental data and mathematical descriptions can
be matched to a high degree of accuracy.  Nevertheless, there are
significant problems associated with conductance-based neuron models
that can make them frustrating to construct and difficult to
interpret.  The most problematic feature is the large number of free
parameters in these model.  For example, the maximal conductances
parameters are typically set by hand until the activity of the model
neuron matches that of its biological counterpart.  Most models are
quite sensitive to the precise values of the maximal conductances, and
experimental measurements cannot specify them accurately enough to
avoid this fine-tuning procedure.  The dynamics of these models is
typically complex, and matching a given type of activity can be a
frustrating exercise in parameter fitting.  Given the difficulty that
modelers have in finding the right parameters to generate a particular
pattern of activity, it is natural to ask how neurons are able to
develop and maintain the conductances they need to function properly.

Ion channels in neurons are continually being synthesized, transported
to various parts of the cell and inserted into the membrane while old
channels are removed and disassembled.  While in the membrane,
channels are subject to a variety of modulatory influences.  It seems
likely that the electrical activity of a neuron plays a role in
regulating all of these processes, so that the neuron can actively
maintain the ion channels it needs to function properly.  Thinking
along these lines, LeMasson, Marder and Abbott (1993) proposed a model
in which the maximal conductances of membrane currents are regulated
by the activity of the neuron through the intracellular \Ca\ 
concentration. In this model, the maximal conductances are not fixed
parameters as in conventional models, but instead can change over time
governed by the intracellular \Ca\ level, [\Ca ].  Because \Ca\ enters
a neuron through voltage-dependent channels, the intracellular \Ca\ 
concentration is a good indicator of electrical activity.  Thus, the
regulation of maximal conductances by [\Ca ] provides a feedback loop
linking the electrical activity of a neuron to the conductances that
produce it.  Bell (1992) and Stemmler and Koch (1999) have proposed
related models and provide information-theoretic interpretations of
activity-dependent conductance regulation.

\begin{figure}
\epsfig{figure=figs/actDepMod.eps, width=4.0in}
\caption{\small{The effect of changing $E_{\Kmy}$ on a model neuron with activity-dependent
    conductances.  Initially the model was in a bursting state.  At
    the time indicated by the triangle, a rise in extracellular \K\ 
    was simulated by changing the value of $E_{\Kmy}$.  This caused
    the neuron to switch its firing pattern.  After a readjustment
    period, bursting was restored. (Adapted from LeMasson {\it et al},
    1993)}}
\label{B2:actDepMod} 
\end{figure}
Neuron models with activity-dependent regulation of conductances can
self-assemble, adjusting their maximal conductances until a fixed
pattern of activity has been established.  Once assembled, these
models are extremely stable and robust.  Figure \ref{B2:actDepMod}
shows a model that is initially in a stable configuration in which it
fired action potentials in bursts. At the time indicated by the
triangle in figure \ref{B2:actDepMod}, an increase in the amount of
extracellular \K\ was simulated by changing the \K\ reversal potential
$E_{\Kmy}$.  This caused the neuron to shift into a fast firing mode,
rather than bursting.  Due to the increased activity, the
intracellular \Ca\ concentration increased and the maximal
conductances were modified by the model.  The result of this
modification was a return to a bursting mode of activity as seen in
the bottom trace of figure \ref{B2:actDepMod}.  Thus, models with
activity-dependent conductances are far more robust than those with
fixed conductances, and they may crudely approximate the types of
adaptations real neurons are capable of making in order to develop and
maintain their electrical properties.

\subsection*{Integrate-and-Fire Models}

When the membrane potential of a neuron reaches a threshold value,
normally between -55 and -50 mV, it will typically fire an action
potential. The action potential follows a rapid stereotyped
trajectory, and then the membrane potential returns to a value
hyperpolarized relative to the action potential threshold.  As we have
seen, the mechanisms by which voltage-dependent \K\ and \Na\ 
conductances produce action potentials are well-understood and can be
modeled quite accurately.  On the other hand, neuron models can be
simplified, and simulations can be accelerated dramatically, if the
biophysical mechanisms responsible for the action potentials are not
explicitly included in the model.  This is the approach used in
integrate-and-fire models.  Integrate-and-fire models generate action
potentials by imposing a fixed rule rather than by computing the
membrane potential trajectory on the basis of modeled conductances.
These models simply stipulate that an action potential occurs whenever
the membrane potential of the model neuron reaches a threshold value
$V_{\thmy}$.  After that action potential, the potential is reset to a
value $V_{\reset}$.

The basic integrate-and-fire model was proposed by Lapicque in 1907,
long before the mechanisms that generate the action potential were
understood.  Despite their age and simplicity, integrate-and-fire
models are still an extremely useful description of neuronal activity.
By avoiding a biophysical description of the action potential,
integrate-and-fire models are left with the simpler task of modeling
only subthreshold membrane potential dynamics.  This can be done with
various levels of rigor.  In the simplest version of these models, all
active membrane conductances are ignored, and the entire membrane
conductance is modeled as a single passive leakage term, $i_{\m} =
\overline g_{\Lmy}(V - E_{\Lmy})$.  For small fluctuations about the
resting membrane potential, neuronal conductances are approximately
constant, and the passive integrate-and-fire model assumes that this
constancy holds over the entire subthreshold range.  For some neurons
this is a reasonable approximation, and for others it is not.  With
these approximations, the model neuron behaves like an electronic
circuit consisting of a resistor and a capacitor in parallel (figure
\ref{B2:ifcircuit}A).
\begin{figure}[hhtt!] 
\epsfig{figure=figs/ifcircuit.eps, width=4.5in, clip}
\caption{\small{The passive integrate-and-fire model.  A) The equivalent resistor-capacitor
    circuit.  B) A passive integrate-and-fire model driven by a
    time-varying current.  The upper trace is the membrane potential
    and the bottom trace the driving current.  The action potentials
    in this figure are simply pasted onto the membrane potential
    trajectory when it reaches the threshold value.  The parameters of
    the model are $E_{\Lmy} = -65$ mV, $V_{\thmy} = -50$ mV,
    $V_{\reset} = -70$ mV, $\tau_{\m} = 10$ ms, and $R_{\m} = 10$
    M$\Omega$.}}
\label{B2:ifcircuit}
\end{figure}

For the time being, we will ignore synaptic inputs by setting $i_{\s}$
to zero and concentrate on the effect of electrode currents.  In this
case, the membrane potential in the passive integrate-and-fire model
is described by the basic equation of a single compartment model,
\ref{B2:onec}, with a single passive membrane conductance,
\begin{equation}
c_{\m}\frac{dV}{dt} = -\overline g_{\Lmy}(V - E_{\Lmy}) + \frac{I_{\emy}}{A}\, .
\end{equation}
It is convenient to multiply this equation by the specific membrane
resistance $r_{\m}$, which is given by $r_{\m} = 1/\overline g_{\Lmy}$
because the model has only a leakage conductance. This cancels the
factor of $\overline g_{\Lmy}$ on the right side of the equation and
leaves a factor $c_{\m}r_{\m} = \tau_{\m}$ on the left side, where
$\tau_{\m}$ is the membrane time constant of the neuron.  The
electrode current ends up being multiplied by $r_{\m}/A$ which is the
total membrane resistance $R_{\m}$.  Thus, the basic equation of the
passive integrate-and-fire models is \begin{equation}
\tau_{\m}\frac{dV}{dt} = E_{\Lmy} - V + R_{\m}I_{\emy} \, .
\label{B2:iandf}
\end{equation}
Equation \ref{B2:iandf} is augmented by the rule that when $V$ reaches
the threshold value $V_{\thmy}$, an action potential is fired and the
potential is reset to $V_{\reset}$.  Equation \ref{B2:iandf} indicates
that when $I_{\emy} = 0$, the membrane potential relaxes exponentially
with time constant $\tau_{\m}$ to $V = E_{\Lmy}$.  $E_{\Lmy}$ is thus
the resting potential of the model cell.  Figure \ref{B2:ifcircuit}B
shows an example of an integrate-and-fire neuron driven by a time
varying external current.

The firing rate of an integrate-and-fire model in response to a
constant injected current can be computed analytically (methods for
integrating equation \ref{B2:iandf} both analytically and numerically
are discussed in Appendix A).  When $I_{\emy}$ is independent of time,
the subthreshold potential $V(t)$ determined by equation
\ref{B2:iandf} is
\begin{equation}
V(t) = E_{\Lmy} + R_{\m}I_{\emy} + (V(0) - E_{\Lmy} - R_{\m}I_{\emy})\exp(-t/\tau_{\m})
\label{B2:soln}
\end{equation}
where $V(0)$ is the value of $V$ at time $t=0$.  This solution can be
checked simply by substituting it into equation \ref{B2:iandf}, but it
is valid for the full integrate-and-fire model only as long as $V$
stays below the threshold.  Suppose that at $t=0$, the neuron has just
fired an action potential and is thus at the reset potential, so that
$V(0) = V_{\reset}$.  The next action potential will occur when the
membrane potential reaches the threshold, that is, at a time $t = T$
when \begin{equation} V(T) = V_{\thmy} = E_{\Lmy} + R_{\m}I_{\emy} +
(V_{\reset} - E_{\Lmy} - R_{\m}I_{\emy})\exp(-T/\tau_{\m}) \, .
\end{equation}
By solving this for $T$, the time of the next action potential, we can determine the
interspike interval, or equivalently, its inverse, the firing rate of the neuron for
constant $I_{\emy}$,
\begin{equation}
r = \frac{1}{T} = \left[\tau_{\m}\ln\left(\frac{R_{\m}I_{\emy} + E_{\Lmy} - V_{\reset}
}{R_{\m}I_{\emy}  + E_{\Lmy} - V_{\thmy}}\right)\right]^{-1} \, .
\end{equation}
This expression is valid if $R_{\m}I_{\emy} > V_{\thmy} - E_{\Lmy}$,
otherwise the firing rate is zero.  For large values of $I_{\emy}$, we
can use the linear approximation of the logarithm ($\ln (1+z) \approx
z$ for small $z$) to show that $r$ grows linearly with $I_{\emy}$,
$r\approx R_{\m}I_{\emy}/(\tau_{\m}(V_{\thmy} - V_{\reset}))$.
 
\begin{figure}[hhtt!]
\epsfig{figure=figs/ifrate.eps, width=4.5in}
\caption{\small{A) Comparison of firing rates as a function of injected current for an
    integrate-and-fire model and a cortical neuron measure {\it in
      vivo}.  The line gives the firing rate of a model neuron with
    $\tau_{\m} = 30$ ms, $V_{\reset} = V_{\rest} = -65$ mV, $V_{\thmy}
    = -50$ mV and $R_{\m} = 90$ M$\Omega$.  The data points are from a
    pyramidal cell in the primary visual cortex of a cat.  The filled
    circles show the inverse of the interspike interval for the first
    two spikes fired, while the open circles show the steady-state
    firing rate after spike-rate adaptation. B) A recording of the
    firing of a cortical neuron under constant current injection
    showing spike-rate adaptation.  C) Membrane voltage trajectory and
    spikes for an integrate-and-fire model with an added current with
    $r_{\m}\Delta g_{\sra}$ = 0.06, $\tau_{\sra}$ = 100 ms, and
    $E_{\Kmy}$ = -70 mV\@.  Other parameters are as in figure
    \ref{B2:ifcircuit}. (Data in A from Ahmed {\it et
al}, 1997, B from McCormick, 1990)}}
%Shephard p. 326
\label{B2:ifrate}
\end{figure}
Figure \ref{B2:ifrate}A compares the firing rate as a function of
$I_{\emy}$, using appropriate parameter values, with data from current
injection into a cortical neuron {\it in vivo}.  The firing rate of
the cortical neuron in figure \ref{B2:ifrate}A has been defined as the
inverse of the interval between pairs of spikes.  The rates determined
in this way using the first two spikes fired by the neuron in response
to the injected current (filled circles in figure \ref{B2:ifrate}A)
agree fairly well with the results of the integrate-and-fire model
with the parameters given in the figure caption.  However, the real
neuron exhibits spike-rate adaptation, the interspike intervals
lengthen over time when a constant current is injected into the cell
(figure \ref{B2:ifrate}B) before settling to a steady-state value.
The steady-state firing rate in figure \ref{B2:ifrate}A (open circles)
could also be fit by an integrate-and-fire model, but not using the
same parameters as were used to fit the initial spikes.  Not all
neurons show spike-rate adaptation, but consideration of this
phenomenon allows us to show how the integrate-and-fire model can be
modified to incorporate more complex dynamics.

\subsubsection*{Spike-Rate Adaptation and Refractoriness}

The passive integrate-and-fire model that we have described is based
on two approximations, a highly simplified description of the action
potential and a linear approximation for the total membrane current.
If details of the action potential generation process are not
important for a particular modeling goal, the first approximation can
be retained while the membrane current is modeled in as much detail as
is necessary.  We will illustrate this process by developing a
heuristic description of spike-rate adaptation using a model
conductance that has characteristics similar to measured neuronal
conductances known to play important roles in producing this effect.

We model spike-rate adaptation by including an addition current in the
model consider thus far,
\begin{equation}
\tau_{\m}\frac{dV}{dt} = E_{\Lmy} - V - r_{\m}g_{\sra}(V - E_K)+ R_{\m}I_{\emy} \, .
\end{equation}
The spike-rate adaptation conductance $g_{\sra}$ has been modeled as a
\K\ conductance so, when activated, it will hyperpolarize the neuron,
slowing any spiking that may be occurring.  We assume that this
conductance relaxes to zero exponentially with time constant
$\tau_{\sra}$ through the equation
\begin{equation}
\tau_{\sra}\frac{dg_{\sra}}{dt} = - g_{\sra} \, .
\end{equation}
Whenever the neuron fires a spike, $g_{\sra}$ is increased by an
amount $\Delta g_{\sra}$, that is, $g_{\sra}\rightarrow g_{\sra} +
\Delta g_{\sra}$. During repetitive firing, the current builds up in a
sequence of steps causing the firing rate to adapt.  Figures
\ref{B2:ifrate}B and \ref{B2:ifrate}C compare the output of the model
with the adapting firing pattern of a cortical neuron.

After a neuron fires an action potential, there is a short interval,
called the absolute refractory period, during which it cannot fire a
second spike.  This is followed by a longer period, the relative
refractory period, when action potentials are more difficult to evoke.
Refractory effects are not included in the basic integrate-and-fire
neuron, but they can be incorporated by adding a conductance like the
spike-rate adaptation conductance, but with a faster recovery time and
a larger conductance increment following an action potential.  With a
large increment, the current can essentially clamp the neuron to
$E_{\Kmy}$ following a spike and temporarily prevent further firing.
As this conductance relaxes back to zero, firing will be possible, but
initially less likely, due to its effects.  When recovery is
completed, normal firing can resume.  An alternative scheme that is
sometimes used to model refractory effects is to raise the threshold
for action potential generation following a spike and then to allow it
to relax back to its normal value.

\subsubsection*{Adding Synaptic Conductances}

Up to now, we have ignored synaptic conductances in the models we have
been discussing.  Synaptic inputs are incorporated into a
conductance-based model by including synaptic conductances (as
described in chapter 5) in the synaptic current appearing in equation
\ref{B2:onec}.  Equivalently, for an integrate-and-fire model a
synaptic conductance term is added to equation \ref{B2:iandf},
\begin{equation}
\tau_{\m}\frac{dV}{dt} = E_{\Lmy} - V - r_{\m}\overline g_{\s}s(V - E_{\s}) + R_{\m}I_{\emy} 
\, .
\label{B2:iandfs}
\end{equation}
The synaptic conductance is written as the product of a fixed maximal
conductance $\overline g_{\s}$and a dynamic gating variable $s$.  The
synaptic current is multiplied by $r_{\m}$ in equation \ref{B2:iandfs}
because equation \ref{B2:iandf} was multiplied by this factor.  To
model synaptic transmission, $s$ changes whenever the presynaptic
neuron fires an action potential.  Various ways of describing and
computing $s$ are discussed in chapter 5\@.

\begin{figure}
\epsfig{figure=figs/IandF2.eps, width=4.5in}
\caption{\small{Two synaptically coupled integrate-and-fire neurons. A) Excitatory
    synapses ($E_{\s}$ = 0 mV) produce an alternating, out-of-phase
    pattern of firing.  B) Inhibitory synapses ($E_{\s}$ = -80 mV)
    produce synchronous firing.  Both model neurons have $V_{\rest}$ =
    -70 mV, $V_{\thmy}$ = -54 mV, $V_{\reset}$ = -80 mV,
    $r_{\m}\overline g_{\s}$ = 0.05, $R_{\m}I_{\emy}$ = 25 mV, and
    $\tau_{\s}$ = 10 ms.}}
\label{B2:iandf2}
\end{figure} 
Figures \ref{B2:iandf2}A and \ref{B2:iandf2}B show examples of two
integrate-and-fire neurons connected by identical excitatory or
inhibitory synapses.  The synaptic conductances in this example are
described by the $\alpha$ function model discussed in chapter 5\@.
This means that the synaptic conductance a time $t$ after the
occurrence of a presynaptic action potential is given by $s =
(t\tau_s)\exp(-t/\tau_s)$ (see chapter 5 for methods of implementing
this model).  The figure shows an interesting effect first noted by
Rinzel and Wang (1992) and further analyzed by Van Vreeswijk,
Ermentrout, and Abbott (1994).  When the synaptic time constant is
sufficiently long ($\tau_{\s}$ = 10 ms in this example), excitatory
connections produce a state in which the two neurons fire alternately,
out of phase with each other.  Inhibitory synapses produces
synchronous firing.  This is a fairly general, though somewhat
counterintuitive, phenomenon.  Inhibitory connections are often more
effective than excitatory connections at synchronizing neuronal firing
(see exercise xx for further analysis).

Synapses have multiple effects on their postsynaptic targets.  To
identify these effects it is useful to express the membrane potential
in terms of a deviation $v$ from its resting value by writing $V =
E_{\Lmy} + v$.  Expressed as an equation for $v$, equation
\ref{B2:iandfs} then becomes
\begin{equation}
\frac{\tau_{\m}}{1 + r_{\m}\overline g_{\s}s}\frac{dv}{dt} = - v + \frac{r_{\m}\overline
g_{\s}s(E_{\s} - E_{\Lmy}) + R_{\m}I_{\emy}}{1 + r_{\m}\overline g_{\s}s} \, .
\, .
\end{equation}
For constant $I_{\emy}$, $v$ approaches the stead-state value
\begin{equation}
v_\infty = \frac{r_{\m}\overline g_{\s}s(E_{\s} - E_{\Lmy}) + R_{\m}I_{\emy}}{1 +
r_{\m}\overline g_{\s}s}
\label{B2:vinfty}
\end{equation}
exponentially with time constant $\tau_{\m}/(1 + r_{\m}\overline
g_{\s}s)$.  This result illustrates three different, but related,
effects of the synapse.  The term proportional to the synaptic
conductance in the numerator of equation \ref{B2:vinfty} represents a
source of current that is added to the external current.  The term in
the denominator of equation \ref{B2:vinfty} shows a divisive effect of
the synapse.  This term also appears in the effective time constant
$\tau_{\m}/(1 + r_{\m}\overline g_{\s}s)$, which is lengthenned by the
synaptic conductance.  Some inhibitory synapses have $E_{\s} \approx
E_{\Lmy}$ so their dominant effect in equation \ref{B2:vinfty} is
divisive.  Such synapses are called shunting, and they provide a
possible biophysical basis for the divisive normalization of visual
responses discussed in chapter 2\@.  Shunting synapses have also been
used in a number of models when a divisive rather than a subtractive
effect of inhibition is needed.

\subsubsection*{Regular and Irregular Firing Modes}

Integrate-and-fire models are useful for studying how neurons
integrate large numbers of synaptic inputs and how networks of neurons
interact.  Softky and Koch (1992 \& 1994) raised an interesting issue
concerning the variability of the response of an integrate-and-fire
neuron receiving a large number of synaptic inputs.  This led to the
realization that neurons can respond to multiple synaptic inputs in
two different modes of operation (Troyer and Miller, 1997; Bugmann
{\it et al}, 1997) depending on the balance that exists between
excitatory and inhibitory inputs (Shadlen and Newsome, 1994).

\begin{figure}
  \epsfig{figure=figs/ifmode.eps, width=4.5in}
\caption{\small{The regular and irregular firing modes of an integrate-and-fire model
    neuron.  A) The regular firing mode.  Upper panel: The membrane
    potential of the model neuron when the spike generation mechanism
    is turned off.  The average membrane potential is above the
    spiking threshold (dashed line).  Lower panel: When the spike
    generation mechanism is turned on, it produces a regular spiking
    pattern. B) The irregular firing mode.  Upper panel: The membrane
    potential of the model neuron when the spike generation mechanism
    is turned off.  The average membrane potential is below the
    spiking threshold (dashed line).  Lower panel: When the spike
    generation mechanism is turned on, it produces an irregular
    spiking pattern.  In order to keep the firing rates comparable in
    these two examples, the value of the reset voltage is higher and
    the membrane time constant is smaller in B than in A\@.}}
\label{B2:modes}
\end{figure} 
The two modes of operation are illustrated in figure \ref{B2:modes},
which shows membrane potentials of an integrate-and-fire model neuron
responding to 1000 excitatory and 200 inhibitory inputs.  Each input
consists of an independent Poisson spike train driving a synaptic
conductance.  The upper panels of figure \ref{B2:modes} show the
membrane potential with the action potential generation mechanism of
the model turned off, and figures \ref{B2:modes}A and \ref{B2:modes}B
illustrate the two different modes of operation.  In figure
\ref{B2:modes}A, the effect of the excitatory inputs is strong enough,
relative to that of the inhibitory inputs, to make the average
membrane potential, when action potential generation is blocked, more
depolarized than the spiking threshold of the model (the dashed line
in the figure).  When the action potential mechanism is turned on
(lower panel of figure \ref{B2:modes}A), this produces a fairly
regular pattern of action potentials.

The irregularity of a spike train can be quantified using the
coefficient of variation ($C_V$), the ratio of the standard deviation
to the mean of the interspike intervals (see chapter 1).  For the
Poisson inputs being used in this example, $C_V$ = 1\@, while for the
spike train in the lower panel of figure \ref{B2:modes}A, $C_V$ =
0.16.  Thus, the output spike train is much more regular than the
input trains.  This is not surprising, because the model neuron
effectively averages its many synaptic inputs.  In the regular firing
mode, the total synaptic input attempts to charge the neuron above the
threshold, but every time the potential reaches the threshold it gets
reset and starts charging again.  In this mode of operation, the
timing of the action potentials is determined primarily by the
charging rate of the cell, which is controlled by its membrane time
constant.

Figure \ref{B2:modes}B shows the other mode of operation that produces
an irregular firing pattern.  In the irregular firing mode, the
average membrane potential is more hyperpolarized than the threshold
for action potential generation (upper panel of figure
\ref{B2:modes}B).  Action potentials are only generated when there is
a fluctuation in the total synaptic input strong enough to make the
membrane potential reach the threshold.  This produces an irregular
spike train, such as that seen in the lower panel of figure
\ref{B2:modes}B which has a $C_V$ value of 0.84\@.

The high degree of variability seen in the spiking patterns of {\it in
  vivo} recordings of cortical neurons suggests that they are better
approximated by an integrate-and-fire model operating in an
irregular-firing mode (see chapter 1).  There are advantages to
operating in the irregular-firing mode that may compensate for its
increased variability.  One is that neurons firing in the irregular
mode reflect in their outputs the temporal properties of fluctuations
in their total synaptic input.  In the regular firing mode, the timing
of output spikes is only weakly related to the temporal character of
the input spike trains.  In addition, neurons operating in the
irregular firing mode can respond more quickly to changes in
presynaptic spiking patterns and firing rates than those operating in
the regular firing mode.

\section{Multi-Compartment Models}

The single-compartment models we have discussed to this point cannot
represent the morphological structure of real neurons.  To do this we
need to use multiple compartments.  It is possible to construct models
with many compartments that reflect the morphology of a given neuron
with great accuracy.  Each compartment in such a model may have a
number of different conductances, resulting in a large number of free
parameters, in particular the maximal conductances for the currents
within each compartment.  At present, these are by no means fully
constrained by experimental data.  Nevertheless, more information
concerning the distribution of conductances in real neurons is
becoming available, and with some assumptions about these
distributions, multi-compartment models can provide key insights
concerning the role of morphology and channel distributions in the
generation of neuronal responses.

\begin{figure}[hhtt!]
\epsfig{figure=figs/compcable.eps, width=4.5in, clip}
\caption{\small{A multi-compartment model of a neuron.  The expanded region shows three
    compartments at a branch point where a single cable splits into
    two.  Each compartment has membrane and synaptic conductances, as
    indicated by the equivalent electrical circuit, and the
    compartments are coupled together by resistors.}}
\label{B2:compcable}
\end{figure}
In a multi-compartment model, each compartment has its own membrane
potential $V_\mu$ (where $\mu$ labels compartments), and its own
gating variables that determine the membrane current for compartment
$\mu$, $i_{\m}^{(\mu)}$\@.  Each membrane potential $V_\mu$ satisfies
an equation similar to (\ref{B2:onec}) except that the compartments
couple to their neighbors in the multi-compartment structure.  For a
nonbranching cable, each compartment is coupled to two neighbors, and
the equations for the membrane potentials of the compartments have the
form \begin{eqnarray} c_{\m}\frac{dV_\mu}{dt} &=& -\sum_i g_i^\mu(V -
E_i) + I_{\emy}^\mu/A_\mu \nonumber \\ &&+ g_{\mu, \mu+1}(V_{\mu+1} -
V_\mu) + g_{\mu, \mu-1}(V_{\mu-1} - V_\mu) \, .
\label{B2:multiComp}
\end{eqnarray}
For a compartment at the end of a cable, there is only one neighboring
compartment, and thus only a single term in the bottom line of
equation \ref{B2:multiComp}.  For a compartment where a cable branches
in two, there are three terms corresponding to coupling of the
branching node to the first compartment in each daughter branch.

The constant $g_{\mu, \mu'}$ that determines the resistive coupling
between neighboring compartments $\mu$ and $\mu'$ are determined by
computing the current that flows between two neighboring compartments
due to Ohm's law.  For simplicity, we begin by computing the coupling
between two compartment that have the same length $L$ and radius $a$.
Using the results of chapter 5, the resistance between two such
compartments, measured from their centers, is the intracellular
resistivity, $r_{\Lmy}$ times the distance between the compartment
centers divided by the cross-sectional area, $r_{\Lmy}L/\pi a^2$\@.
The total current flowing between compartment $\mu+1$ and compartment
$\mu$ is then $\pi a^2 (V_{\mu+1}- V_\mu)/r_{\Lmy}L$\@.  Recall that
equation \ref{B2:onec} for the potential within a single compartment
refers to currents per unit area of membrane.  Thus, in the equation
for compartment $\mu$ of a multi-compartment model we must divide the
total current from compartment $\mu'$ by the surface area of
compartment $\mu$, $2\pi aL$\@.  Thus, we find that $g_{\mu, \mu'} =
a/(2r_{\Lmy}L^2)$.

The value of $g_{\mu, \mu'}$ is given by a more complex expression if
the two neighboring compartments have different lengths or radii.
This can occur when a tapering cable is approximated by a sequence of
cylindrical compartments, or at a branch point where a single
compartment connects with two other compartments as in figure
\ref{B2:compcable}.  In either case, suppose that compartment $\mu$
has length $L_\mu$ and radius $a_\mu$ and compartment $\mu'$ has
length $L_{\mu'}$ and radius $a_{\mu'}$.  The resistance between these
two compartments is the sum of the two resistances from the middle of
each compartment to the junction between them, $r_{\Lmy}L_\mu/(2\pi
a_\mu^2) + r_{\Lmy}L_{\mu'}/(2\pi a_{\mu'}^2)$.  To compute
$g_{\mu,\mu'}$ we divide this by the total surface area of compartment
$\mu$, $2\pi a_\mu L_\mu$, which gives
\begin{equation}
g_{\mu,\mu'} = \frac{a_\mu a_{\mu'}^2}{r_{\Lmy}L_\mu (L_\mu a_{\mu'}^2 +
L_{\mu'}a_\mu^2)} \, .
\end{equation}

Equations such as \ref{B2:multiComp} for all of the compartments of a
model neuron determine the membrane potential throughout the neuron
with a spatial resolution given by the compartment size.  An efficient
method for integrating equation \ref{B2:multiComp} is discussed in
Appendix C\@.  Using this scheme, models can be numerically integrated
extremely efficiently, even those involving large numbers of
compartments.  Such integration schemes are built into neuron
simulation software packages such as Neuron and Genesis.

\subsection*{Action Potential Propagation Along an Unmyelinated Axon}

As an example of multi-compartment modeling, we simulate the
propagation of an action potential along an unmyelinated axon using
the currents measured in the squid giant axon by Hodgkin and Huxley.
In the Hodgkin-Huxley model, $i_{\m}$ is the sum of three currents, a
leakage current, a delayed-rectifier \K\ current, and a transient \Na\ 
current, and these depend on three dynamic variables $m$, $h$ and $n$.
Specifically,
\begin{equation} 
i_{\m} ={\overline g}_{\Lmy}(V-E_{\Lmy})+{\overline g}_{\Kmy}n^4(V-E_{\Kmy})+
{\overline g}_{\Namy}m^3h(V-E_{\Namy})
\end{equation} 
where ${\overline g}_L = 3$ $\mu$S/mm$^2$, ${\overline g}_K = 360$
$\mu$S/mm$^2$, ${\overline g}_{Na} = 1200$ $\mu$S/mm$^2$, $E_{\Lmy}$ =
-54.402 mV, $E_{\Kmy}$ = -77 mV and $E_{\Namy}$ = 50 mV\@.  The
dynamic equations for the variables $m$, $h$, and $n$ are given in
chapter 5\@.  Figure \ref{B2:apaxon} shows an action potential
propagating along an axon with each compartment modeled in this way.
The action potential extends over more than 1 mm of axon and it
travels at a speed of about 2 mm/ms or 2 m/s.
\begin{figure}[hhtt!]
\epsfig{figure=figs/apaxon.eps, width=4.5in, clip}
\caption{\small{Propagation of an action potential along a multi-compartment model axon. 
    The upper panel shows the multi-compartment representation of the
    axon with 100 compartments.  The axon segment shown is 4 mm long
    and has a radius of 1 $\mu$m.  An external current sufficient to
    initiate action potentials is injected at the point marked
    $I_{\emy}$.  The panel beneath this shows the membrane potential
    as a function of position along the axon, at a given instant of
    time.  The spatial position in this panel is aligned with the axon
    depicted above it.  The action potential is moving to the right.
    The bottom two panels show the membrane potential as a function of
    time at the two locations denoted by the arrows and symbols $V_1$
    and $V_2$ in the upper panel.}}
\label{B2:apaxon}
\end{figure}

\subsection*{Propagation Along a Myelinated Axon}

Many axons in vertebrates are covered with an insulating sheath of
myelin except at gaps, called the nodes of Ranvier, where
voltage-dependent ion channels (primarily \Na\ channels) are
concentrated (see figure \ref{B2:myaxon}A).  Action potentials
propagate passively down the myelin-covered sections of the axon and
are actively regenerated at the nodes of Ranvier.  The myelin sheath
consists of many layers of membrane wrapped around the axon.  This
greatly increases the membrane resistance and greatly decreases the
capacitance of the covered segment of the axon.  Figure
\ref{B2:myaxon} shows the equivalent circuit for a multi-compartment
model of a myelinated axon.
\begin{figure}[hhtt!]
\epsfig{figure=figs/meylin.eps, width=4.5in}
\caption{\small{A myelinated axon.  A) The equivalent circuit for a multi-compartment
    representation of a myelinated axon.  The myelinated segments are
    represented by a membrane capacitance, a passive membrane
    resistance, and a longitudinal resistance.  The nodes of Ranver
    contain voltage-dependent conductances (primarily \Na\ with some
    \K ) as well.  B) A cross-section of a myelinated axon consisting
    of a central axon core of radius $a_1$ and a myelin sheath making
    a total outside radius of $a_2$.}}
\label{B2:myaxon} 
\end{figure}

We can compute the membrane resistance of a myelin covered axon by
treating the myelin sheath as an extremely thick cell membrane.
Consider the geometry shown in the cross-sectional diagram of figure
\ref{B2:myaxon}B\@.  The myelin sheath extends from the radius $a_1$
of the axon core to the outer radius $a_2$.  If the thickness of a
single layer of cell membrane is $d_{\m}$, the specific membrane
resistance per unit membrane thickness is $r_{\m}/d_{\m}$, where
$r_{\m}$ is the specific membrane resistance introduced in chapter
5\@.  The membrane resistance of a cylindrical shell of membrane with
thickness $\Delta a$, radius $a$, and length $\Delta x$ is the
specific membrane resistance per unit thickness times the thickness of
the shell divided by the area of the shell, which is $r_{\m}\Delta
a/(d_{\m}2\pi a\Delta x)$.  The total membrane resistance of a length
$\Delta x$ of myelinated axon is obtained by taking the limit $\Delta
a\rightarrow 0$ and integrating over the thickness of the sheath,
\begin{equation} 
R_{\m} = \frac{r_{\m}}{d_{\m}2\pi \Delta x}\int_{a_1}^{a_2}\frac{da}{a} = 
\frac{r_{\m}\ln(a_2/a_1)}{d_{\m}2\pi \Delta x} \, . 
\end{equation}
A similar calculation of the membrane capacitance for this length of
axon gives the result $C_{\m} = c_{\m}d_{\m}2\pi \Delta
x/\ln(a_2/a_1)$, so we see that the membrane resistance and membrane
capacitance are scaled by inverse factors.  This means that the
membrane time constant $\tau_{\m} = R_{\m}C_{\m} = r_{\m}c_{\m}$ is
not affected by the myelination.  The longitudinal resistance of this
length of axon, $R_{\Lmy} = r_{\Lmy}\Delta x/(\pi a_1^2)$, is also
unaffected by the myelination, because longitudinal current flows only
inside the central core.
 
The calculations done in the previous paragraph can be used to
determine the optimal thickness of myelination for an axon of fixed
outside diameter.  Specifically, we can determine the core diameter
$a_1$ that gives the maximum action potential propagation speed if an
axon is required to fit within a space of radius $a_2$.  As shown in
chapter 5, the effective speed of propagation down a passive cable is
$2\lambda/\tau_{\m}$, where $\lambda$ is the length constant of the
cable.  Because myelination has no effect on the membrane time
constant $\tau_{\m}$, we must choose $a_1$ to make the length constant
$\lambda$ as big as possible to maximize this propagation speed.  In
chapter 5, we showed that the length constant for an unmyelinated axon
was given by $\lambda^2 = ar_{\m}/(2r_{\Lmy})$.  A review of this
derivation shows that the generalization to the myelinated case is
\begin{equation} \lambda^2 = \frac{(\Delta x)^2 R_{\m}}{R_{\Lmy}} =
  \frac{r_{\m}a_1^2\ln(a_2/a_1)}{2d_{\m}r_{\Lmy}} \, .
\end{equation}
To determine the inner radius that maximizes this length constant, we
set the derivative of this expression with respect to $a_1$ to zero.
This gives the result $a_1 = a_2\exp(-1/2)$ or $a_1 = 0.6a_2$.  An
inner core fraction of 0.6 is typical for myelinated axons, indicating
that the thickness of myelin is indeed optimal.
  
\section{Firing-Rate Models}

The models we have discussed thus far produce spike sequences
deterministically in response to injected current or synaptic input.
Sequences of action potentials produced by these models can only be
compared with those of real neurons if the total synaptic input to the
real neuron is known, and this is rarely the case.  If the synaptic
input is not known, it may only be possible to perform statistical
comparisons of a model with the data.  For example, a model may fail
to predict exact spike sequences, although it adequately predicts
firing rates.  In this case, a simpler model that is couched in terms
of firing rates rather than spike sequences may provide a sufficient
description.  Such models were discussed in chapters 1 and 2, where we
also discussed methods of generating spike trains stochastically from
them.

Single neuron models are the basic building blocks for constructing
models of neural circuits and networks.  In a network model, we are
primarily concerned with faithfully reproducing the dynamics of the
circuit being modeled.  Network dynamics is affected both by the
response properties of individual neurons and by the synaptic
connections between neurons.  Thus, in addition to reproducing
individual neuronal responses, a single neuron model used within a
network must generate realistic synaptic inputs to other neurons.
Here models that produce sequences of spikes have a distinct advantage
because the spike sequences can be fed into a model of synaptic
conductances to produce realistic synaptic input.  Nevertheless,
reduced models that deal exclusively with firing rates, rather than
generating spike sequences, are used extensively in network modeling
because of their ease of analysis and because, in some regimes at
least, they are adequately similar to their spiking counterparts.  For
example, firing-rate models form the basis for the discussions of
network models, adaptation, and learning in the following chapters.
In this section, we introduce firing-rate models of neural networks
and discuss the assumptions upon which they are based. These models
make most sense in the case that there are many, interconnected,
neurons. We therefore treat this case formally. In the next chapter we
consider some computational properties of the resulting dynamical
system. 

Consider a network of $N$ neurons that will be described by firing
rates $r_a$ with $a=1, 2, \ldots , N$.  If neuron $b$ fires at time
zero, we write the post-synaptic potential observed at the cell body
of neuron $a$ (the somatic potential) at time $t$ as $M_{ab}K(t)$.
Writing it like this implies that we are ignoring any
voltage-dependence in the conductances and the dendrites.  $M_{ab}$
describes the amplitude and sign of the post-synaptic potential
generated in neuron $a$ by neuron $b$, and $K(t)$ describes its time
course.  To avoid an ambiguity in the definition of $M_{ab}$, we
normalize $K(t)$ by requiring its integral over all positive times to
be one.  For simplicity at this point, we use the same function $K(t)$
to describe all synapses.

Assuming linearity, the total somatic potential of neuron $a$ at time
$t$ due to a sequence of presynaptic spikes at synaptic terminal $b$
occurring at times $t_i$ is the sum
\begin{equation}
M_{ab}\sum_{t_i<t}K(t - t_i)  =
M_{ab}\intl{0}{\infty}{\tau}K(\tau)\rho_b(t-\tau) \, .
\end{equation}
On the right side of this equation, we have used the neural response
function introduced in chapter 1, $\rho_b(t) = \sum_i\delta(t - t_i)$,
to describe the sequence of spikes fired by neuron $b$.  The equality
follows from integrating over the sum of $\delta$ functions in the
definition of $\rho_b(t)$.  The total synaptic contribution to the
potential at the cell body of neuron $a$, $V_a(t)$ is obtained by
summing the contributions from all of its synaptic inputs,
\begin{eqnarray}
V_a(t) & = & \sum_{b=1}^N
M_{ab}\intl{0}{\infty}{\tau}K(\tau)\rho_b(t-\tau) \,  \label{B2:totsyn}
\\[5pt]
& = & \intl{0}{\infty}{\tau}K(\tau) \sum_{b=1}^N
M_{ab} \rho_b(t-\tau) \ .\label{eq:totsynp} 
\end{eqnarray}
In a spiking model such as the integrate-and-fire model we considered
above, the synaptic contribution to this potential should be augmented
by a term to take account of the reset of the cell to $V_\reset$
following the emission of a spike. The spike response model (see
Gerstner, 1998), which is a generalisation of the integrate-and-fire
model, takes explicit account of these extra terms in order to study
aspects of the very short-term behavior of networks of neurons.
Standard firing rate models ignore it, and instead make two critical
approximations. One is to replace $\rho_b(t)$, which is a sum of a set of
delta functions, with its trial-based average
\begin{equation}
\rho_b(t) \approx \langle \rho_b(t) \rangle = p_b(t)
\end{equation}
where $p_b(t)$ is the instantaneous firing rate of neuron $b$.
The second approximation is that the firing rate of neuron $b$,
in the form of the
instantaneous spike probability, $p_b(t)$, is related to the membrane
potential by the instantaneous equation:
\begin{equation}
p_b(t) = F(V_b(t))
\label{eq:pform}
\end{equation}
where function $F$ will be specified later.  Some conditions under
which these approximations are both close and far are known. Very
crudely, they are safe away from either very fast firing rates or very
short-term behavior of the cells (such as synchronization of
individual spikes).  They rely on a form of mean-field approximation
used in statistical physics, that the variability across a large
number of inputs effectively reproduces the trial-to-trial variability
for any single input, so the sum over synaptic inputs effectively
averages out the trial-to-trial variability on any one of them.

The two forms of $V_a(t)$ in
equations~\ref{B2:totsyn}~and~\ref{eq:totsynp} lead to two slightly
different dynamical equations for the evolution of the rates, which
are based on different assumptions about the longer of two timescales,
for the {\it synapse\/} or for the {\it membrane\/} (Ermentrout,
1998). The form in equation~\ref{B2:totsyn} is best interpreted in the
case that the dynamics of the membrane are fast compared with those of
the synapse, or the cell sits so close to its spiking threshold that
little integration by the membrane is required. In this case, the
membrane principally acts to add the contributions from the synapses,
which themselves are synaptically filtered versions of the presynaptic
activity. Write
\begin{equation}
U_b(t) = \intl{0}{\infty}{\tau}K(\tau)p_b(t-\tau)
= \intl{-\infty}{t}{\tau}K(t-\tau)p_b(\tau) \ ,
\label{eq:Uform}
\end{equation}
for the integral in equation~\ref{B2:totsyn}.  Note that $U_b(t)$ is
an appropriate candidate for $r_b(t)$, which is a temporally {\it
  smoothed\/} version of the firing rate of neuron $b$, and so we use
$r_b(t)$ in its place.  Applying
equations~\ref{eq:pform}~and~\ref{eq:Uform} to neuron $a$, we get
\begin{equation}
r_a(t) = \intl{-\infty}{t}{\tau}K(t-\tau)p_a(\tau)
= \intl{-\infty}{t}{\tau}K(t-\tau)F\left(\sum_{b=1}^NM_{ab} r_b(\tau)\right) 
\label{eq:Ufulleq}
\end{equation}
If a synapse can be treated as a simple low pass filter with a
time constant $\tau_s$, {\it ie\/}
\begin{equation}
K(\tau) = \frac{1}{\tau_s} e^{-\frac{\tau}{\tau_s}}
\end{equation}
then we can differentiate equation~\ref{eq:Ufulleq} to get
\begin{equation}
\tau_s \frac{d}{dt} r_a(t) = - r_a(t) + F\left(\sum_{b=1}^N M_{ab}
r_b(t)\right) 
\label{B2:frmod}
\end{equation}
If the synaptic contribution to the potential is held constant, $r_a$
approaches the value $F(\sum M_{ab}r_b)$, so $F$ is the steady-state
firing rate as a function of synaptic current.  It is often modeled as
a threshold linear function $F(\sum M_{ab}r_b) = [\sum M_{ab}r_b -
\Theta]_+$, where $\Theta$ is the threshold voltage and, as in
previous chapters, the notation $[\:]_+$ denotes rectification.
Alternatively, $F$ is often modeled as a saturating function such as a
sigmoidal function.  Sometimes an external input $h_a$ to neuron $a$
is included by writing $F(h_a + \sum M_{ab}r_b)$ in equation
\ref{B2:frmod}.

The form of the integral equation~\ref{eq:totsynp} suggests an
alternative interpretation in the case that the synapses are fast
compared with the membrane. Then, the synaptic input
$\sum_{b=1}^N M_{ab} \rho_b(t)$ is filtered by the membrane
$K(\tau)=\frac{1}{\tau_m} e^{-\frac{\tau}{\tau_m}}$ with the membrane
time constant $\tau_m$, leading, using equation~\ref{eq:pform}, to
\begin{equation}
V_a(t) = \intl{-\infty}{t}{\tau} K(t-\tau) \sum_{b=1}^N M_{ab}
F\left(V_b(\tau)\right) \ .
\end{equation}
Differentiating this gives
\begin{equation}
\tau_m \frac{d}{dt} V_a(t) = -V_a(t) + \sum_{b=1}^N M_{ab}
F\left(V_b(t)\right) 
\label{eq:ptlform}
\end{equation}
The firing rate of cell $a$ is $p_a(t)=F(V_a(t))$. Since
equation~\ref{eq:ptlform} is essentially a low pass filter, again we
can equate this with a smoothed firing rate $r_a(t)$.
Equation~\ref{eq:ptlform} is sometimes called the membrane potential
form of the firing rate equation, whereas equation~\ref{B2:frmod} is
called the activity form. Multiplying equation~\ref{B2:frmod} by
matrix $M$ and writing $\tilde{V}_a(t)=\sum_{b=1}^N M_{ab} r_b(t)$, we
get
\begin{equation}
\tau_s \frac{d}{dt} \tilde{V}_a(t) = -\tilde{V}_a(t) + \sum_{b=1}^N
M_{ab} F\left(\tilde{V}_b(t)\right)
\label{eq:reduceptl}
\end{equation}
just as in equation~\ref{eq:ptlform} (though with a different time
constant). The dynamical behaviors of the two equations are therefore
very closely related.

\subsection*{Paired Excitatory and Inhibitory Neurons}

Some models use a single {\it pair\/} of an excitatory and an
inhibitory cell as their basic computational unit. This can take
better account of the neural separation between excitation and
inhibition, even though such strict pairing is not actually present.
Following Wilson \& Cowan (1972), we analyse the range of basic
dynamical behaviors of such pairs using the firing rate formulation.
We consider the case of a single pair, with the firing rate of the
excitatory cell being $r_{\Emy}$ and the inhibitory cell being
$r_{\Imy}$; although this might more reasonably come from a population
of pairs, with identical firing rates.

In this example, the excitatory and inhibitory synapses are allowed to
have different time constants, $\tau_{\Emy}$ and $\tau_{\Imy}$, and
the excitatory and inhibitory firing rates are defined using
exponential kernels with these time constants.  The recurrent synaptic
strength from the excitatory neuron to itself is denoted by
$M_{\Emy\Emy}$, from the inhibitory neuron to itself by
$M_{\Imy\Imy}$, from the inhibitory to the excitatory neuron by
$M_{\Emy\Imy}$, and from the excitatory to the inhibitory neuron by
$M_{\Imy\Emy}$.  Using threshold linear response functions, and
equation~\ref{B2:frmod}, the firing rates are determined by
\begin{equation}
\tau_{\Emy}\frac{dr_{\Emy}}{dt} = -r_{\Emy} +\left[M_{\Emy\Emy}r_{\Emy} +
M_{\Emy\Imy}r_{\Imy} - \Theta_{\Emy}\right]_+
\label{B2:pop1}
\end{equation}
and
\begin{equation}
\tau_{\Imy}\frac{dr_{\Imy}}{dt} = -r_{\Imy} + \left[M_{\Imy\Imy}r_{\Imy} +
M_{\Imy\Emy}r_{\Emy} - \Theta_{\Imy}\right]_+ \, .
\label{B2:pop2}
\end{equation}
In the example we consider, we set $M_{\Emy\Emy} = 1.25$,
$M_{\Imy\Emy} = 1$, $M_{\Imy\Imy} = M_{\Emy\Imy} = -1$, $\Theta_{\Emy}
= -10$ Hz, $\Theta_{\Imy} = 10$ Hz, $\tau_{\Emy} = 10$ ms, and we vary
the value of $\tau_{\Imy}$.  The negative value of $\Theta_{\Emy}$
means that this parameter serves as a source of background activity
(activity in the absence of excitatory input) rather than as a
threshold.

\subsection*{Phase-Plane Methods and Stability Analysis}

The model of the pair of interacting excitatory and inhibitory cells
given by equations \ref{B2:pop1} and \ref{B2:pop2} provides an
opportunity to illustrate some of the techniques used to study the
dynamics of nonlinear systems.  This model exhibits both static
(constant $r_{\Emy}$ and $r_{\Imy}$) and oscillatory activity
depending on the values of its parameters.  Stability analysis can be
used to determine the parameter values where transitions between these
two types of activity take place.

The firing rates $r_{\Emy}(t)$ and $r_{\Imy}(t)$ arising from
equations \ref{B2:pop1} and \ref{B2:pop2} can be displayed by plotting
them as functions of time, as in figures \ref{B2:EIss}A and
\ref{B2:EIosc}A\@.  Another useful way of depicting these results,
illustrated in figures \ref{B2:EIss}B and \ref{B2:EIosc}B, is to plot
pairs of points $(r_{\Imy}(t), r_{\Emy}(t))$ for all relevant $t$
values.  As the firing rates change, these points trace out a curve or
trajectory in the $r_{\Emy}$-$r_{\Imy}$ plane, which is called the
phase plane of the model.  Phase-plane plots can be used to give a
geometric picture of the dynamics of a model.

\begin{figure}
\epsfig{figure=figs/lambda.eps, width=4.5in}
\caption{\small{A) Nullclines, flow directions, and fixed point for the firing-rate model of
    interacting excitatory and inhibitory populations.  The two
    straight lines are the nullclines along which $dr_{\Emy}/dt = 0$
    or $dr_{\Imy}/dt = 0$.  The filled circle is the fixed point of
    the model.  The horizontal and vertical arrows indicate the
    directions that $r_{\Emy}$ (horizontal arrows) and $r_{\Imy}$
    (vertical arrows) flow in different regions of the phase plane
    relative to the nullclines.  B) Real (upper panel) and imaginary
    (lower panel) parts of the eigenvalue determining the stability of
    the fixed point.  To the left of the point where the imaginary
    part of the eigenvalue goes to zero, both eigenvalues are real.
    The imaginary part has been divided by $2\pi$ to give the
    frequency of oscillations near the fixed point.}}
\label{B2:lambda}
\end{figure} 
The values of $r_{\Emy}$ and $r_{\Imy}$ for which the right sides of
either equation \ref{B2:pop1} or equation \ref{B2:pop2} vanish are of
particular interest in phase-plane analysis.  These sets of values
form two curves in the phase plane known as nullclines.  The
nullclines for equations \ref{B2:pop1} and \ref{B2:pop2} are the
straight lines drawn in figure \ref{B2:lambda}A\@.  The nullclines are
interesting because they divide the phase plane into regions with
opposite flow patterns.  This is because the right side of equation
\ref{B2:pop1} or equation \ref{B2:pop2} is positive on one side of its
nullcline and negative on the other.  The right sides of these
equations determine the sign and magnitude of the rate of change of
the corresponding firing rate, and thus set the direction of flow in
the phase plane.  These flow directions are denoted by horizontal and
vertical arrows in figure \ref{B2:lambda}A\@.  Thus, above the
nullcline along which $dr_{\Emy}/dt = 0$, $dr_{\Emy}/dt < 0$, and
below it $dr_{\Emy}/dt > 0$.  Similarly, $dr_{\Imy}/dt > 0$ to the
right of the nullcline where $dr_{\Imy}/dt = 0$, and $dr_{\Imy}/dt <
0$ to the left of it.

A fixed point is a pair of firing rates that make $dr_{\Emy}/dt =
dr_{\Imy}/dt =0$.  Since both derivatives vanish, these values will
remain fixed in time.  Because a fixed point requires both derivatives
to vanish, it can only occur at an intersection of nullclines.  The
model we are considering has a single fixed point denoted by the
filled circle in figure \ref{B2:lambda}A\@.  A fixed point provides a
potential static configuration for the system, but it is critically
important whether the fixed point is stable or unstable.  If a fixed
point is stable, initial values of $r_{\Emy}$ and $r_{\Imy}$ near the
fixed point will be drawn toward it over time.  If the fixed point is
unstable, nearby configurations are pushed away from the fixed point,
and the system will only remain at the fixed point indefinitely if the
rates are set initially to the fixed-point values with infinite
precision.

\begin{figure}
\epsfig{figure=figs/EIss.eps, width=4.5in}
\caption{\small{Activity of the excitatory-inhibitory firing-rate model when the
    fixed point is stable.  A) The excitatory and inhibitory firing
    rates settle to the fixed point over time.  B) The phase-plane
    trajectory is a counter-clockwise spiral collapsing to the fixed
    point. For this example, $\tau_{\Imy} = 30$ ms, and all other
    parameters are as given in the text.}}
\label{B2:EIss}
\end{figure} 
Linear stability analysis (see chapter 12) can be used to determine
whether the fixed point is stable or unstable.  This analysis starts
by considering the first derivatives of the right sides of equations
\ref{B2:pop1} and \ref{B2:pop2} with respect to $r_{\Emy}$ and
$r_{\Imy}$ evaluated at the values of $r_{\Emy}$ and $r_{\Imy}$
corresponding to the fixed point.  The four combinations of
derivatives computed in this way can be arranged into a matrix
\begin{equation}
\left(
\begin{array}{cc}
(M_{\Emy\Emy} - 1)/\tau_{\Emy} 
& M_{\Emy\Imy}/\tau_{\Emy}\\
M_{\Imy\Emy}/\tau_{\Imy} 
& (M_{\Imy\Imy}-1)/\tau_{\Imy}
\end{array}
\right) \, .
\label{eq:blockm}
\end{equation}
As discussed in chapter 12, the stability of the fixed point is
determined by the eigenvalues of this matrix.  The computation of
these eigenvalues is discussed in chapter 12\@.  They are given by
\begin{equation}
\lambda = \frac{1}{2}\left(\frac{M_{\Emy\Emy} - 1}{\tau_{\Emy}}+
\frac{M_{\Imy\Imy}-1}{\tau_{\Imy}} \pm \sqrt{\left(\frac{M_{\Emy\Emy} - 1}{\tau_{\Emy}}-
\frac{M_{\Imy\Imy}-1}{\tau_{\Imy}}\right)^2 +
4\frac{M_{\Emy\Imy}M_{\Imy\Emy}}{\tau_{\Emy}\tau_{\Imy}}}\right)
\label{B2:eigen}
\end{equation}
The stability of the fixed point is determined by the real parts of
these eigenvalues (see chapter 12).  If both real parts are less than
zero the fixed point is stable, while if either is greater than zero
the fixed point is unstable.  If the factor inside the square root in
equation \ref{B2:eigen} is positive, both eigenvalues are real, and
the behavior near the fixed point will be exponential.  This means
that there will either be exponential movement toward the fixed point
(if both eigenvalues are negative) or away from the fixed point if
either eigenvalue is positive.  We focus on the case when the factor
inside the square root is negative, so that the square root is
imaginary and the eigenvalues form a complex conjugate pair. In this
case, the behavior near the fixed point is oscillatory and the
trajectory will either spiral into the fixed point, if the real part
of the eigenvalues is negative, or outward from the fixed point if the
real part of the eigenvalues is positive.  The imaginary part of the
eigenvalue determines the frequency of oscillations near the fixed
point.
\begin{figure}
\epsfig{figure=figs/EIosc.eps, width=4.5in}
\caption{\small{Activity of the excitatory-inhibitory firing-rate model when the
    fixed point is unstable.  A) The excitatory and inhibitory firing
    rates settle into periodic oscillations.  B) The phase-plane
    trajectory is a counter-clockwise spiral that joins the limit
    cycle, which is the closed orbit.  For this example, $\tau_{\Imy}
    = 50$ ms and all other parameters are as given in the text.}}
\label{B2:EIosc}
\end{figure} 
The real and imaginary parts of one of these eigenvalues are plotted
as a function of $\tau_{\Imy}$ in figure \ref{B2:lambda}B\@.  This
figure indicates that the fixed point is stable if $\tau_{\Imy} < 40$
ms and unstable for larger values of $\tau_{\Imy}$.

Figures \ref{B2:EIss} and \ref{B2:EIosc} show examples in which the
fixed point is stable and unstable, respectively. In figure
\ref{B2:EIss}A, the oscillations in $r_{\Emy}$ and $r_{\Imy}$ are
damped, and the firing rates settle down to the stable fixed point
values.  The corresponding phase-plane trajectory is a collapsing
spiral (figure \ref{B2:EIss}B).  In figure \ref{B2:EIosc}A the
oscillations grow, and in figure \ref{B2:EIosc}B the trajectory is a
spiral that expands outward until the system enters a limit cycle.  A
limit cycle is a closed orbit in the phase plane indicating periodic
behavior.  The fixed point is thus unstable in this case, and the
limit cycle is stable.  The transition from a fixed point to limit
cycle behavior as a parameter is changed (in this case $\tau_{\Imy}$)
making the real part of a complex eigenvalue change sign is called a
Hopf bifurcation.  This is one of the ways that a nonlinear system can
make a transition from a stable fixed point to a limit cycle.  In a
Hopf bifurcation, the limit cycle emerges at a finite frequency, which
is similar to the behavior of a type II neuron when it starts firing
action potentials.  Other types of bifurcations produce type I
behavior with oscillations emerging at zero frequency.  One example of
this is a saddle node bifurcation, which occurs when two fixed points,
one stable and one unstable, meet at the same point in phase space.

Firing-rate models greatly simplify neuronal dynamics, but they cannot
account for temporal details of neuronal spiking beyond those
described by a firing rate.  Nevertheless, these models provide a
compact description of at least some aspects of neuronal activity, and
they are used extensively in studies of large neural networks.  We
will consider computational properties of firing-rate models in
considerably more detail in chapter 7\@.


\section{Appendices}

\subsection*{A) Integrating the Membrane Potential} 

We begin by considering the integration of equation \ref{B2:iandf}.  It is convenient to
rewrite this equation in the form
\begin{equation}
\tau_V\frac{dV}{dt} = V_\infty - V  \, .
\label{B2:Veq}
\end{equation}
where $\tau_V = \tau_m$ and $V_\infty = E_{\Lmy} + R_{\m}I_{\emy}$.   When the external
current $I_{\emy}$ is independent of time, the solution of this equation is
\begin{equation}
V(t) = V_\infty + (V(t_0) - V_\infty)\exp(-(t-t_0)/\tau_V)
\label{B2:Vsoln}
\end{equation}
where $t_0$ is any time prior to $t$ and $V(t_0)$ is the value of $V$ at time $t_0$.
Equation \ref{B2:soln} is a special case of this result with $t_0 = 0$.

If $I_{\emy}$ depends on time, the solution \ref{B2:Vsoln} is not valid.  An analytic
solution can still be written down in this case (see exercise xx), but it is not
particularly useful except in special cases.  Over a small enough time period $\Delta t$,
we can approximate $I_{\emy}(t)$ as a constant and use the solution \ref{B2:Vsoln} to step
from a time $t$ to $t + \Delta t$.  This requires replacing the variable $t_0$ in equation
\ref{B2:Vsoln} with $t$ and $t$ with $t + \Delta t$ so that 
\begin{equation}
V(t+\Delta t) = V_\infty + (V(t) - V_\infty)\exp(-\Delta t/\tau_V) \, .
\label{B2:Vstep}
\end{equation}
This equation provides an updating rule for the numerical integration of equation
\ref{B2:Veq}.  Provided that $\Delta t$ is sufficiently small, repeated application of the
update rule \ref{B2:Vstep} provides an accurate way of determining the membrane
potential. Furthermore, this method is stable because, if $\Delta t$ is too large, it will
only move $V$ toward $V_\infty$ and not, for example, make it grow without bound.

The equation for a general single-compartment conductance-based model, equation \ref{B2:onec}
with \ref{B2:memcurr}, can be written in the same form as equation \ref{B2:Veq} with
\begin{equation}
V_\infty = \frac{\sum_i\overline g_im_i^{p_i}h_i^{q_i}E_i + I_{\emy}/A}{\sum_i\overline
g_im_i^{p_i}h_i^{q_i}}
\end{equation} 
and
\begin{equation}
\tau_V = \frac{c_{\m}}{\sum_i\overline g_im_i^{p_i}h_i^{q_i}} \, .
\end{equation}
Note that if $c_m$ is in units of nF/mm$^2$ and the maximal conductances in $\mu$S/mm$^2$,
$\tau_V$ comes out in ms units.  Similarly, if the reversal potentials are given in units of
mV, $I_{\emy}$ is in nA, and $A$ is in mm$^2$, $V_{\infty}$ will be in mV units.  

If we take the time interval $\Delta t$ to be small enough so that the gating variables
can be approximated as constant during this period, the membrane potential can again be
integrated over one time step using equation \ref{B2:Vstep}.  Of course, the gating
variables are not fixed, so once $V$ has been updated by this rule, the gating variables
must be updated as well.

\subsection*{B) Integrating the Gating Variables} 

All the gating variables in a conductance-based model satisfy equations of the same form,
\begin{equation} 
\tau_z\frac{dz}{dt} = z_\infty - z
\end{equation}
where we use $z$ to denote a generic variable.  Note that this equation has the same
form as equation \ref{B2:Veq}, and it can be integrated in exactly the same way.  We assume
that $\Delta t$ is sufficiently small that $V$ does not change appreciably over
this time interval (and similarly [\Ca ] is approximated as constant over this interval if
any of the conductances are \Ca -dependent).  Then, $\tau_z$ and $z_\infty$, which are
functions of $V$ (and possibly [\Ca ]) can be treated as constants over this period and $z$
can be updated by a rule identical to
\ref{B2:Vstep},
\begin{equation} 
z(t+\Delta t) = z_\infty + (z(t) - z_\infty)\exp(-\Delta t/\tau_z).
\label{B2:zstep}
\end{equation}
An efficient integration scheme for conductance-based models is to alternate using rule
(\ref{B2:Vstep}) to update the membrane potential and rule (\ref{B2:zstep}) to update all
the gating variables.  It is important to alternate the updating of $V$ with that of the
gating variables, rather than doing them all simultaneously, as this keeps the method
accurate to second order in $\Delta t$\@.  If \Ca -dependent conductances are included,
the intracellular \Ca\ concentration should be computed simultaneously with the
membrane potential.

\subsection*{C) Integrating Multi-Compartment Models}

To solve the compartmental equations numerically, we must integrate one set of
equations for each compartment.  The gating variables can be handled as in the single
compartment case using equation (\ref{B2:zstep}). Integrating the membrane potentials for
the different compartments is more complex because they are coupled to each other.

The equation for the membrane potential within compartment $\mu$, equation
\ref{B2:multiComp}, can be written in the form 
\begin{equation} 
\frac{dV_\mu}{dt} = A_\mu V_{\mu-1} + B_\mu V_\mu + C_\mu V_{\mu+1} + D_\mu 
\label{B2:comp}
\end{equation}
where $A_\mu = c_m^{-1}g_{\mu, \mu-1}$, $B_\mu = -c_m^{-1}(\sum_i g_i^\mu + g_{\mu, \mu+1} +
g_{\mu, \mu-1})$, $C_{\mu} = c_m^{-1}g_{\mu, \mu+1}$, and $D_{\mu} = c_m^{-1}(\sum_i g_i^\mu
E_i + I_{\emy}^\mu/A_\mu)$.  Note that the gating variables and other parameters have been
absorbed into the values of $A_\mu$, $B_\mu$, $C_\mu$, and $D_\mu$ in this equation. 
Equation \ref{B2:comp}, with $\mu$ running over all of the compartments of the model,
generates a set of coupled differential equations.  Because of the coupling between
compartments, we cannot use the rule \ref{B2:Vstep} to integrate these equations. Instead,
we present another method that shares some of the positive features of rule \ref{B2:Vstep}.

Two of the most important features of an integration method are accuracy and stability. 
Accuracy refers to how closely numerical finite-difference methods reproduce the exact
solution of a differential equation as a function of the integration step size $\Delta t$. 
Stability refers to what happens when $\Delta t$ is chosen to be excessively large and the
method starts to become inaccurate.  A stable integration method will degrade
smoothly as $\Delta t$ is increased, producing results of steadily decreasing accuracy.  An
unstable method, on the other hand, will, at some point, display a sudden transition and
generate wildly inaccurate results.  For the present purposes, we want a method that is
stable, even if this costs us somewhat in accuracy.  

Defining
\begin{equation}  
V_\mu(t+\Delta t) = V_\mu(t) + \Delta V_\mu \, ,
\end{equation}
the finite difference form of equation \ref{B2:comp} gives the update rule
\begin{equation} 
\Delta V_\mu  = \left(A_\mu V_{\mu-1}(t) + B_\mu V_\mu(t) + C_\mu V_{\mu+1}(t) +
D_\mu\right)\Delta t
\label{B2:euler}
\end{equation}
which is how $\Delta V_\mu$ is computed using the so-called Euler method.  This method is
both inaccurate and unstable.  The stability of the method can be improved dramatically by
evaluating the membrane potentials on the right side of equation \ref{B2:euler} not at time
$t$, but at a later time $t + \rho \Delta t$, so that
\begin{equation}
\Delta V_\mu = [A_\mu V_{\mu-1}(t+\rho \Delta t) + B_\mu V_\mu(t+\rho \Delta t)
 + C_\mu V_{\mu+1}(t+\rho \Delta t) + D_\mu]\Delta t \, .
\label{B2:deltaV}
\end{equation}
Two such methods are predominantly used, the reverse Euler method for which $\rho=1$ and the
Crank-Nicholson method with $\rho=0.5$\@.  The reverse Euler method is the more stable of
the two and the Crank-Nicholson is the more accurate.  In either case, $\Delta
V_\mu$ is determined from equation \ref{B2:deltaV}.  These methods are called implicit
because equation \ref{B2:deltaV} must be solved to determine $\Delta V_\mu$.  To do this, we
write $V_\mu(t+\rho \Delta t) \approx V_\mu(t) + \rho \Delta V_\mu$ and likewise for
$V_{\mu\pm 1}$.  Substituting this into equation \ref{B2:deltaV} gives
\begin{equation}
\Delta V_\mu = a_\mu\Delta V_{\mu-1} + b_\mu\Delta V_\mu + c_\mu\Delta V_{\mu+1} + d_\mu
\label{B2:implicit}
\end{equation}
where $a_\mu = A_\mu\rho\Delta t$, $b_\mu = B_\mu\rho\Delta t$, $c_\mu = C_\mu\rho\Delta
t$ and $d_\mu = [D_\mu + A_\mu V_{\mu-1}(t) + B_\mu V_\mu(t) + C_\mu V_{\mu+1}(t)]\Delta t$.

Equation \ref{B2:implicit} for all $\mu$ values provides a set of coupled linear equations
for the needed quantities $\Delta V_\mu$\@.  An efficient method exists for solving these
equations (Hines 1984, Tuckwell 1988).  We illustrate the method for a single, unbranched
cable and leave the extension to a branched cable as an exercise for the reader (exercise
xx).  The cable we consider begins with the compartment $a=1$, so $a_1 = 0$ and it ends at
compartment $N$, so $c_N = 0$\@.  The method consists of solving (\ref{B2:implicit}) for
$\Delta V_\mu$ in terms of  $\Delta V_{\mu+1}$ sequentially starting at one end of the cable
and proceeding to the other end.  For example, if we start the procedure at compartment one,
$\Delta V_1$ can be expressed as \begin{equation} \Delta V_1 = \frac{c_1\Delta V_2 +
d_1}{1-b_1} . \end{equation} Substituting this into the equation (\ref{B2:implicit}) for
$\Delta V_2$ gives \begin{equation}
\Delta V_2 = \tilde{b}_2\Delta V_2 + c_2\Delta V_3 + \tilde{d}_2
\end{equation}
where $\tilde{b}_2 = b_2 + a_2c_1/(1-b_1)$ and $\tilde{d}_2 = d_2 + a_2d_1/(1-b_1)$\@.  We
now repeat the procedure going down the cable. At each stage, we solve for $\Delta V_{\mu-1}$
in terms of $\Delta V_\mu$ finding
\begin{equation}
\Delta V_{\mu-1} = \frac{c_{\mu-1}\Delta V_\mu + \tilde{d}_{\mu-1}}{1-\tilde{b}_{\mu-1}}  .
\label{B2:veqn}
\end{equation}
and define
\begin{equation}
\tilde{b}_{\mu+1} = b_{\mu+1} + \frac{a_{\mu+1}c_\mu}{1-\tilde{b}_\mu}
\label{B2:btilde}
\end{equation}
and
\begin{equation}
\tilde{d}_{\mu+1} = d_{\mu+1} + \frac{a_{\mu+1}\tilde{d}_\mu}{1-\tilde{b}_\mu}  .
\label{B2:dtilde}
\end{equation}
Finally, when we get to the end of the cable we can solve for
\begin{equation}
\Delta V_N = \frac{\tilde{d}_N}{1-\tilde{b}_N}
\label{B2:N}
\end{equation}
because $c_N = 0$\@.  

The procedure for finding the $\Delta V_\mu$ is the following.   Define $\tilde {b_1} = b_1$
and $\tilde{d}_1 = d_1$ and iterate equations \ref{B2:btilde} and \ref{B2:dtilde} down
the length of the cable to define all the $\tilde{b}$'s and $\tilde{d}$'s.  Then, solve for
$\Delta V_N$ from equation \ref{B2:N} and iterate back up the cable solving for the $\Delta
V$'s using \ref{B2:veqn}\@.  This process takes only $2N$ steps.  As the reader can show
in exercise xx, a similar procedure can be used for a branching cable.  The equations are
solved starting at the ends of both branches and moving in toward the branching node, then
continuing on as for an unbranched cable, and finally reversing direction and completing the
solution moving in the opposite direction along the cable and i  ts branches. 

\section{Further Reading and References}

\subsubsection{Further Reading}
\begin{list}{}{\usecounter{bean}
\setlength{\leftmargin}{0.2in}
\setlength{\listparindent}{-0.2in}}
\item \mbox{ }\vspace{-0.8cm}

Hodgkin, AL \& Huxley, AF (1952) A quantitative description of membrane
current and its application to conduction and excitation in nerve.\ {\it Journal of Physiology
(London)} {\bf 117}:500-544.

Johnston, D \& Wu, SM (1995) {\it Foundations of Cellular Neurophysiology}.\ Cambridge, MA:MIT Press.

Koch, C (1998) {\it Biophysics of Computation: Information Processing in Single Neurons}.\ New York:Oxford
University Press.

Koch, C \& Segev, I, editors (1998) {\it Methods in Neuronal Modeling: From Synapses to Networks}.\
Cambridge, MA:MIT Press.

McCormick, DA (1990) Membrane properties and neurotransmitter actions.\ In GM Shepherd,
editor, {\it The Synaptic Organization of the Brain}.\ New York:Oxford University Press.

Tuckwell, HC (1988) {\it Introduction to Theoretical Neurobiology}.\ Cambridge, UK:Cambridge
University Press.

\end{list}

\subsubsection{References}
\begin{list}{}{\usecounter{bean}
\setlength{\leftmargin}{0.2in}
\setlength{\listparindent}{-0.2in}}
\item \mbox{ }\vspace{-0.8cm}

Bell, A (1992) Self-Organization in real neurons: anti-Hebb in `channel space'? In
JE Moody and SJ Hanson, editors, {\it Neural Information Processing Systems 4}.\
San Mateo CA:Morgan Kaufman, pp. 59-66.

Bressloff, PC \& Coombes, S (1999) Dynamics of strongly coupled
spiking neurons. Submitted to {\it Neural Computation.\/}

Bugmann, G, Christodoulou, C \& Taylor, JG (1997) Role
of temporal integration and  fluctuation detection in the highly irregular firing of a leaky
integrator neuron model with partial reset.\ {\it Neural Computation} {\bf 9}:985-1000.

Connor, JA \& Stevens CF (1971) Prediction of repetitive firing behavior from voltage clamp
data on an isolated neuron soma.\ {\it Journal of Physiology} {\bf 213}:31-54.

Connor, JA, Walter, D \& McKown, R (1977) Neural repetitive firing: modifications of the
hodgkin-Huxley axon suggested by experimental results from crustacean axons.\ {\it
Biophysical Journal} {\bf 18}:81-102.

Ermentrout, B (1998) Neural networks as spatio-temporal
pattern-forming systems. {\it Reports on Progress in Physics\/} {\bf
  64}:353-430. 

FitzHugh, R (1961) Impulses and physiological states in models of nerve membrane.\ {\it
Biophysical Journal} {\bf 1}:445-466.

Freeman, WJ \& Schneider, W (1982) Changes in spatial patterns of rabbit olfactory
EEG with conditioning to odors.\ {\it Psychophysiology\/} {\bf 19}:44-56.

Gerstner, W (1998) Spiking neurons.  In W Maass \& CM Bishop, editors,
{\it Pulsed Neural Networks\/} Cambridge, MA: MIT press: 3-54.

Haberly, LB (1990) Olfactory cortex.\ In GM Shepherd, editor, {\it The
  Synaptic Organization of the Brain}.\ New York:Oxford University
Press.

Hines, M (1984) Efficient computation of brancehd nerve equations.\ 
{\it International Journal of Biomedical Computation} {\bf 15}:69-76.

Huguenard, JR \& McCormick, DA (1992) Simulation of the currents
involved in rhythmic oscillations in thalamic relay neurons.\ {\it
  Journal of Neurophysiology} {\bf 68}:1373-1383.

LeMasson, G, Marder, E \& Abbott, LF (1993) Activity-Dependent Regulation of
Conductances in Model Neurons.\ {\it Science} {\bf 259}:1915-1917.

Li, Z (1995). Modeling the sensory computations of the olfactory
bulb. In JL van Hemmen, E Domany \& K Schulten, editors, {\it 
Models of Neural Networks, Volume 2.\/} New York: Springer Verlag.
 
Li, Z (1990) A model of olfactory adaptation and sensitivity
enhancement in the olfactory bulb.\ {\it Biological Cybernetics} {\bf
  62}:349-361.

Li, Z \& Hopfield, J (1989) Modeling the olfactory bulb and its neural oscillatory
processings.\ {\it Biological Cybernetics} {\bf 61}:379-392. 

Morris, C \& Lecar, H (1981) Voltage oscillations in the barnacle giant muscle fiber.\ {\it
Biophysical Journal} {\bf 35}:193-213.

Nagumo, JS, Arimoto, S \& Yoshizawa, S (1962) {\it Proceedings of the IRE} {\bf
50}:2061-2074.

Shadlen, MN \& Newsome, WT (1994) Noise, neural codes and cortical organization.\ {\it
Current Opinion in Neurobiology} {\bf 4}:569-579.

Softky, WR \& Koch, C (1992) Cortical cells should spike regularly but do
not.\ {\it Neural Computation} {\bf 4}:643-646.  

Softky, WR \& Koch, C (1994) The highly irregular firing of cortical
cells is inconsistent with temporal integration of random EPSPs.\ {\it
Journal of Neuroscience} {\bf 13}:334-350.

Stemmler, M \& Koch, C (1999) Information maximization in single neurons.\ (unpublished).

Troyer, TW \& Miller, KD (1997) Physiological gain leads to high ISI variability in a simple 
model of a cortical regular spiking cell.\ {\it Neural Computation} {\bf 9}:971-983.

Van Vreeswijk, C, Abbott, LF \& Ermentrout, GB (1994) When Inhibition Not Excitation
Synchronizes Neuronal Firing.\ {\it Journal of Computational Neuroscience} {\bf 1}:313-321.

Wang, X.-J. (1994) Multiple dynamical modes of thalamic relay neurons: Rhythmic bursting
and intermittent phase-locking.\ {\it Neuroscience} {\bf 59}:21-31.

Wang, X.-J. \& Rinzel, J (1992) Alternating and synchronous rhythms in reciprocally
inhibitory model neurons.\ {\it Neural Computation} {\bf 4}:84-97.

Turrigiano, G, LeMasson, G \& Marder, E (1995) Selective regulation of current
densities  underlies spontaneous changes in the activity of cultured neurons.\ {\it Journal
of Neuroscience} {\bf 15}:3640-3652.

Wilson, HR \& Cowan, JD (1972) Excitatory and inhibitory interactions in localized
populations of model neurons.\ {\it Biophysical Journal} {\bf 12}:1-24.

\end{list}

\section{Exercises}
\begin{enumerate}

\item Show that the Hodgkin-Huxley model produces a type II firing-rate curve in response to
constant current injection.

\item Demonstrate postinhibitory rebound in the Hodgkin-Huxley model.

\item Construct a Morris-Lecar model neuron (Morris and Lecar, 1981). This is probably
the simplest conductance-based model with interesting dynamics.  Instead of modeling the
fast sodium spikes of an action potential, it models slower calcium spikes that can produce
bursting behavior in neurons.  The model has just two active currents, an instantaneous
voltage-dependent \Ca\ current and a persistent potassium current, described by a single
dynamical gating variable $n$,
\begin{equation} 
i_{\m} ={\overline g}_L(V-E_{\Lmy})+ {\overline g}_{Ca}m_\infty(V)(V-E_{Ca}) + 
{\overline g}_Kn(V-E_{\Kmy})
\end{equation} 
with ${\overline g}_L$ = 0.5 mS/cm$^2$, ${\overline g}_{Ca}$ = 1 mS/cm$^2$ and
${\overline g}_{K}$ = 2 mS/cm$^2$, $E_{\Lmy}$ = -50 mV, $E_{Ca}$ = 100 mV and $E_{\Kmy}$ =
-70 mV\@.
The function $m_\infty$ is given by
\begin{equation}
m_\infty = \frac{1}{1 + \exp[-.133(V+1)]}
\end{equation}
and the gating variable $n$ is given by
\begin{equation}
\tau_n\frac{dn}{dt} = n_\infty - n 
\end{equation}
with
\begin{equation}
\tau_n = \frac{3}{\cosh[.0345(V -10)]} \;\;\;\;\;
n_\infty = \frac{1}{1+\exp [-.138(V-10)]}
\end{equation}
Determine the firing rate as a function of injected current and plot the membrane
potential both as a function of time and in a phase-plane trajectory along with the $n$
variable.  In the phase plane, plot the nullclines for the $V$ and $n$ equations.

\item Construct a reciprocal two-neuron oscillator model (due to F. Skinner and E. Marder).
For both cells, the membrane current is 
\begin{equation} 
i_{\m} ={\overline g}_L(V-E_{\Lmy})+{\overline g}_Hz(V-E_H)
\end{equation} 
where ${\overline g}_L$ = 0.0111 mS/cm$^2$, ${\overline g}_H$ = 0.5 mS/cm$^2$, $E_{\Lmy}$
= -60 mV, and $E_H$ = -10 mV\@.  The H-current is a hyperpolarization activated
current determined by
\begin{equation}
\tau_z(V)\frac{dz}{dt} = z_\infty(V) - z
\end{equation}
with
\begin{equation}
r_\infty(V) = \frac{1}{1 + \exp\left[(V + 50)/7\right]} \:\:\:\:\mbox{and}\:\:\:\:
\tau_r(V) = \frac{3030.3}{1 + \exp\left[(- V - 110)/13\right]} \, .
\end{equation}
The synaptic connection between these neurons is modeled as a graded function of the
presynaptic membrane potential. The synaptic current for a postsynaptic neuron with membrane
potential $V_{post}$ is
\begin{equation}
I_{syn} ={\overline g}_ss(V_{post}-E_{\s})
\end{equation}
where $E_{\s}$ = -80 mV and ${\overline g}_s$ = 0.5 mS/cm$^2$.  The variable $s$ is
determined by the membrane potential of the presynaptic neuron, $V_{pre}$, by
\begin{equation}
\tau_{\s}(1 - s_\infty(V_{pre}))\frac{ds}{dt} = s_\infty(V_{pre}) - s
\end{equation}
with $\tau_{\s} = 50$ ms and 
\begin{equation}
s_\infty(V_{pre}) = \left[\tanh\left((V_{pre} + 57)/10)\right)\right]_+
\end{equation}
where $[z]_+ = z$ if $z\geq 0$ and $[z]_+ = 0$ if $z < 0$.

\item The FitzHugh-Nagumo equations (FitzHugh, 1961; Nagumo {\it et al}, 1962) are given by
\begin{equation}
\frac{dv}{dt} = v(1-v^2) - u + I_{\emy} \;\;\;\;\; \mbox{and}\;\;\;\;\;
\frac{du}{dt}   = \epsilon (v - 0.5u)
\end{equation}
Draw the nullclines for these equations for $I_{\emy}=0$ and $I_{\emy}=-1$.  These are the
lines in the $x$-$y$ plane where the right sides of these two equations are zero.  In which
case or cases do you think the model will produce oscillations?  Next simulate the model to
see what happens when these equations are integrated over time.  Determine what happens for
$I_{\emy}=0$ with $\epsilon$ = 0.3, 0.1, and 1 and for $I_{\emy}=-1$ with $\epsilon = 0.3$.

\item Write down the analytic solution of equation \ref{B2:Veq} when $I_{\emy}(t)$ is an
arbitrary function of time.  The solution will involve integrals that cannot be performed
unless  $I_{\emy}(t)$ is specified.

\item Construct the model of two, coupled integrate-and-fire model neurons of figure
\ref{B2:iandf2}.  Show how the pattern of firing for the two neurons depends on the
strength, type (excitatory or inhibitory), and time constant of the reciprocal synaptic
connection (see Van Vreeswijk {\it et al}, 1994).

\item  Construct a non-branching axonal cable with conductances in each compartment
described by the Connor-Stevens model.  Plot the action potential propagation velocity
as a function of the axon radius.

\item Simulate the dynamics of the firing-rate model of coupled
  excitatory and inhibitory neural populations given by equations
  \ref{B2:pop1} and \ref{B2:pop2}.  Choose a parameter range different
  from that discussed in the text, and compare the behavior of the
  model to what you expect on the basis of a fixed-point stability
  analysis.

\item Determine the numerical solution for a multi-compartment cable
  with a single branch node analogous to the solution for a
  non-branching cable given in appendix C\@.

\end{enumerate}