 |
INTRODUCTION |
In the absence of stimuli, neurons in many areas of the
mammalian CNS remain active. This activity falls into two
main patterns: low, steady firing in the 1- to 5-Hz range and rhythmic
bursting. The former has been widely observed in sensory cortex, motor
cortex, and spinal cord and has been seen in both awake-behaving and
anesthetized animals (Collins 1987
; Gilbert
1977; Herrero and Headley 1997; Lamour et
al. 1985; Leventhal and Hirsch 1978;
Mednikova and Kopytova 1994; Ochi and Eggermont
1997; Salimi et al. 1994; Szente et al. 1988; Zurita et al. 1994); the latter is
fundamental to central pattern generators, which generate rhythmic
behavior such as respiration, locomotion, and chewing (for reviews see
Marder and Calabrese 1996; Rekling and Feldman
1998). Both patterns may occur in the same neural tissue, with
neuromodulators inducing a reversible switch between steady firing and
bursting (Berkinblit et al. 1978; Kudo and Yamada
1987; Zoungrana et al. 1997).
From a theoretical point of view, it has been relatively easy to
understand how large neuronal networks might generate rhythmic bursting
(Feldman and Cleland 1982; Perkel and Mulloney
1974; Rekling and Feldman 1998; Rybak et
al. 1997; Smith 1997) but difficult to
understand how such networks might generate steady, low firing rates
(Abeles 1991). Several investigators have explored the
latter problem theoretically (Amit and Treves 1989;
Buhmann 1989; Treves and Amit 1989).
Their approach was to use simulated networks and attempt to generate
low firing rates through dynamic interactions between excitatory and
inhibitory neurons. The networks they used were (1) highly
interconnected, (2) isolated from any external sources of
input, and (3) comprised of neurons whose resting membrane potentials were far enough below threshold that several
near-synchronous excitatory postsynaptic potentials (EPSPs) were
necessary to trigger an action potential. Networks with these seemingly
standard properties were unable to produce maintained low firing rates;
the lowest numerically generated mean rates were ~20 Hz,
significantly larger than the background firing rates observed in
biological networks.
These experimental and numerical findings raise three main
questions. First, why were networks with the above "standard"
properties unable to fire at low rates? Did the numerical studies
simply miss a parameter regime that would have generated low firing
rates, or were the properties actually too restrictive? Second, what are the conditions that allow networks to fire robustly at low rates?
And third, what is the relation between steadily firing networks and
rhythmically bursting ones? In particular, is there a single parameter,
or a small set of parameters, that allow networks to switch naturally
from one state to the other?
To answer these questions, we developed a model that allowed us to
explore analytically the intrinsic dynamics of large, isolated networks
in the low firing rate regime. Our main assumptions were conventional:
excitatory input to a neuron increases its firing rate, inhibitory
input decreases it, and neurons exhibit spike-frequency adaptation
(their firing rates decrease during repetitive firing). We found
theoretically, and verified using large, simulated networks of spiking
neurons, that a single parameter plays a dominant role in controlling
intrinsic firing patterns. That parameter is the fraction of
endogenously active cells; i.e., the fraction of cells that fire
without input. Our primary result is that there is a sharp distinction
between networks in which the fraction of endogenously active cells is
zero and those in which it is greater than zero. When the fraction is
zero, a stable, low firing rate equilibrium cannot exist;
networks are either silent or fire at high rates. Only when the
fraction is greater than zero are low average firing rates possible. In
this regime there is a further subdivision into networks that burst and
those that fire steadily, the former having fractions below a threshold
and the latter having fractions above threshold. Connectivity also
plays a role in shaping firing patterns, but to a lesser degree; it
determines whether activity occurs at a constant mean firing rate or
oscillates around that mean.
These theoretical results imply that an isolated network that fires at
low rates must contain endogenously active cells. We tested this
strong, parameter-free prediction experimentally in cultured neuronal
networks and consistently found a large fraction of endogenously active
cells, ~30% on average, in networks that displayed low firing rates.
We also investigated experimentally the transition between steady
firing and bursting as the fraction of endogenously active cells
changed, and we found that networks that fired steadily could be
induced to burst by reducing the fraction of endogenously active cells.
Both of these experimental results are presented in the following paper
(Latham et al. 2000).
The above analysis applies to networks that receive external input as
well as to isolated ones. This is because cells that receive sufficient
external input to make them fire are effectively endogenously active,
in the sense that they can fire without receiving input from other
neurons within the network. These pseudoendogenously active cells
control firing patterns in the same way that truly endogenously active
ones control firing patterns in isolated networks. In particular, for
networks receiving external input, firing patterns follow the same set
of transitions as discussed above: as activity produced by external
input increases there is a transition first from silence to bursting,
and then from bursting to steady firing.
| |
METHODS |
In RESULTS we make a series of predictions about the
behavior of large neuronal networks, then test those predictions by
performing large-scale network simulations. In this section we describe
the single neuron equations and synaptic coupling used in the
simulations, discuss our choice of simulation parameters, and provide a
particular realization of a reduced network model that is used to
illustrate the general principles derived in RESULTS.
Network simulations
Our predictions about the behavior of large neuronal networks
are based on a simplified firing rate model. To verify these predictions, we perform large-scale simulations with synaptically coupled model neurons. Following are the single neuron equations used
in the simulations, the prescription for choosing single neuron
parameters and synaptic coupling, and a complete list of network parameters.
SINGLE NEURON EQUATIONS.
Our simulated network consists of N synaptically coupled
neurons, NE of which are excitatory and
NI of which are inhibitory. The time evolution
equation for the membrane potential of neuron i,
Vi, may be written
|
(1)
|
where Ccell is the cell capacitance,
Ispike consists of the currents responsible for
generating action potentials, IfAHP and
IsAHP are fast and slow afterhyperpolarization
currents, respectively, Isyn represents the
synaptic current, and Ia is a constant
depolarizing current. The fast afterhyperpolarization current enforces
a relative refractory period; the slow afterhyperpolarization current,
which introduces a slowly decaying hyperpolarization each time a neuron emits an action potential, is responsible for spike-frequency adaptation.
Rather than using conductance-based currents for
Ispike, we model this quantity as
|
(2)
|
where 
V 
Vt 
Vr is the nominal gap between resting
(Vr) and threshold (Vt)
voltages and Rcell is the membrane resistance of
the cell (Fig. 1). Typically,
Vr = 65 mV and Vt = 50 mV, producing a nominal gap of 15 mV. In the absence of synaptic
drive and afterhyperpolarization currents, Eq. 1 with
Ispike,i given by Eq. 2 is
identical to the 
-neuron model introduced by Ermentrout and
Kopell (1986) and explored further by Ermentrout (1996). This model describes the behavior of type I neurons,
neurons that can support arbitrarily low-frequency oscillations
(Hodgkin 1948). The advantage of using the -neuron
rather than a conductance-based model is that relatively large time
steps can be used, on the order of 1 ms rather than the 0.1- to 0.2-ms
time steps required to accurately represent the rapidly varying
membrane potential during a true action potential.
View larger version (11K):
[in this window]
[in a new window]
|
Fig. 1.
Current-voltage relation for Ispike, the current
that generates action potentials. Arrows indicate voltage trajectories
in the absence of synaptic input and afterhyperpolarization currents.
Arrows pointing to the right correspond to regions where the voltage
increases in time; arrows pointing to the left correspond to regions
where it decreases. These arrows indicate that
Vr, the resting membrane potential, is a stable
equilibrium and Vt, the threshold, is an
unstable one. Inset: action potential generated by giving
the neuron an initial voltage slightly above threshold.
|
|
If we were to adhere strictly to the -neuron model [which is
related to our model by the change of variables
Vi ~ tan (/2) + constant] an action
potential would occur whenever Vi reaches +
(the spike apex), at which point the voltage would be reset to
(the spike repolarization). For numerical work, however, it is
necessary to use finite spike apex and repolarization voltages. The
apex we use in our simulations, denoted Vapex,
is 20 mV, and the reset value, Vrepol, is 80 mV.
The fast and slow afterhyperpolarization currents may be written
|
(3a)
|
|
(3b)
|
where 
K is the potassium reversal
potential and gK,i and
gK
Ca,i are potassium
conductances. The latter quantities obey the time evolution equations
|
(4a)
|
|
(4b)
|
where gK,
and
gKCa, are the equilibrium value
of the fast and slow potassium conductances,
KCa and
KCa are their time constants,
tiµ is time of the µth spike on neuron
i, and the 
-function, (t tiµ), provides an impulse whenever
t = tiµ. Those impulses
causes discrete increases, gK and
gKCa, in
gK,i and
gKCa,i, respectively.
For simplicity we assume that the equilibrium conductances, gK, and
gKCa,, and the time constants,
KCa and
KCa are independent of voltage.
Then, without loss of generality, we may assume that
gK,i and
gKCa,i are zero at rest,
which implies that gK, = gKCa, = 0.
Combining Eqs. 1-4, the network evolution equations become
|
(5a)
|
|
(5b)
|
|
(5c)
|
where cell RcellCcell is the cell
time constant and quantities with a hat have been multiplied by
Rcell:
Îa,i RcellIa,i,
etc. Note that both Îa,i and
Îsyn,i have units of voltage,
whereas the conductances, 
K and
KCa, and their increments, K and
KCa, are dimensionless.
The normalized synaptic current, Îsyn, is
written
|
(6)
|
where Wij is the strength of the
connection from neuron j to neuron i,
sij(t) is the fraction of channels on
neuron i that open when neuron j emits an action
potential, and j is the reversal potential
associated with the ligand-gated receptor on the cell that is
postsynaptic to neuron j. The weight matrix, Wij, is always nonnegative. This ensures that
excitatory drive (activity from neurons with
j greater than the threshold for the emission
of an action potential) increases the probability that a postsynaptic
neuron fires, whereas inhibitory drive (activity from neurons with
j less than threshold) decreases the probability.
The conductance change on a postsynaptic neuron is mediated by the
fraction of open channels, sij. That fraction
increases instantaneously when a spike occurs at neuron j
and decays exponentially between spikes: for all i
|
(7)
|
The parameter rs determines how many
closed channels open each time neuron j fires. If we had
included adaptation, rs would have been called
the use parameter (Tsodyks and Markram 1997); we adopt
that name here.
Equations 5-7, along with the rule that whenever
Vi reaches Vapex an
action potential is emitted and the voltage is reset to
Vrepol, constitute our complete network
equations. Implementing this model as it stands, however, is costly
numerically. This is because the fraction of open channels,
sij(t), must be integrated separately
for each synapse. These separate integrations can be avoided by noting
that the relevant s-dependent quantity is the sum that
appears on the right hand side of Eq. 6. That sum can be
divided into two terms
|
(8)
|
where
|
(9a)
|
|
(9b)
|
Combining the definitions in Eq. 9 with the time
evolution for the fraction of open channels, Eq. 7, we see
that 
i and 
,i evolve according to
|
(10a)
|
|
(10b)
|
With this formulation there is no need to keep track of the
sij individually; Eq. 8 along with
the time evolution equations for i and
,i, Eq. 10, can
be used to determine the synaptic current, and
sij drops out of the equations.
DISTRIBUTION OF APPLIED CURRENTS.
The applied current, Îa,i,
determines how close a neuron is to the threshold for generating an
action potential, or, if it is above threshold, its endogenous firing
rate. This quantity is chosen probabilistically:
Îa,i has a boxcar distribution,
uniform between 0 and Îmax and 0 outside those two values. We use this distribution because it gives us precise
control over the number of neurons that are endogenously active. The
precise control arises because neuron i is endogenously active if and only if Îa,i > V/4 Îthresh. Consequently, when Îmax < Îthresh no neurons are endogenously active, and when Îmax > Îthresh the fraction of endogenously active cells is (Îmax Îthresh)/Îmax.
NETWORK CONNECTIVITY.
Connectivity is specified by the weight matrix,
Wij. This quantity determines both whether two
neurons are connected and, if they are, the strength of that
connection. Like the applied current, the weight matrix is chosen
probabilistically, and the probability that neuron j
connects to neuron i is denoted Pij. If there is a connection, the strength of that connection is set by
considerations of EPSP and inhibitory postsynaptic potential (IPSP)
sizes. If there is no connection, Wij is set to zero.
We use two models of connectivity, infinite range and local,
when constructing Pij. For infinite range
connectivity, Pij depends only on the types of
neurons i and j
where T specifies type
We do not allow autapses, which means Pii = 0.
The connection probabilities consist of four numbers, corresponding to
the four combinations of pairs of excitatory and inhibitory neurons. We
parameterize these four numbers with the quantities KTj and BTj,
where KTj is the mean number of
postsynaptic neurons a presynaptic neuron connects to and
BTj is the connectivity bias
|
(11a)
|
|
(11b)
|
where, recall, NE and
NI are the number of excitatory and inhibitory
neurons, respectively. Note that BTj > 1 indicates a bias toward inhibitory neurons, whereas
BTj < 1 indicates a bias toward excitatory ones.
Inverting Eq. 11 yields the expressions for the connection
probabilities in terms of the mean number of connections and bias
|
(12a)
|
|
(12b)
|
Equation 12 is sufficient to specify connectivity in
the infinite range regime.
For local connectivity, Pij depends on distance
between neurons. For "distance between neurons" to make sense we
need to assign to each neuron a spatial location. We work in a
two-dimensional space, and specify the location of a neuron in polar
coordinates, (r, ), where r is radius and is azimuthal angle. Neurons are randomly assigned a position according
to the probability distribution P(r, ) where
P(r, )rdrd is the probability
that a neuron falls within the area bounded by [r,
r + dr] and [, + d]. For
P(r, ) we use the azimuthally symmetric
function
![<IT>P</IT>(<IT>r</IT><IT>, &thgr;</IT>)<IT>=</IT><FR><NU><IT>1−tanh </IT>[(<IT>r</IT><SUP><IT>2</IT></SUP><IT>−1</IT>)<IT>/&Dgr;</IT><SUB><IT>r</IT></SUB>]</NU><DE><IT>&pgr;&Dgr;</IT><SUB><IT>r</IT></SUB><IT> ln </IT>[<IT>1+exp</IT>(<IT>2/&Dgr;</IT><SUB><IT>r</IT></SUB>)]</DE></FR>](/content/vol83/issue2/fulltext/808/img027.gif)
|
(13)
|
where r is the width of the transition
region between approximately constant density and near zero density;
i.e., the distribution P(r, ) is approximately
flat for radius r < 1 r and drops rapidly to nearly zero in a
transition region of width 2r.
Using the probability distribution given in Eq. 13 to assign
a position, xi (ri cos i,
ri sin i), to every
neuron i, the distance between two neurons is then given simply by the Euclidean norm, dij |xi xj|. We use a Gaussian profile for local
connectivity, for which the probability of making a connection is
modulated by the factor exp(dij2/2
Tj2)
where Tj is a parameter that determines
the axonal spread of neurons of type j. Thus
where ZTj is chosen so that
the mean connection probability of a neuron placed at the origin
(r = 0) is equal to PTiTj. In
our calculations we use an approximate expression for the normalization, valid in the limit that r 
1 [so that P(r, ) = 1/
when
r 
1 and zero when r > 1]. In this
limit
This expression for ZTj is
valid as long as ZTj 
PTiTj, which
can always be satisfied by taking the limit
Tj 
(in this limit
ZTj 1). We impose the condition ZTj PTiTj in all
our simulations.
We consolidate the infinite range and local connectivity models by
using the local connectivity description and setting
Tj = whenever we want to recapture
the infinite range case.
The strength of a connection, once one is made, is determined by the
size of Wij. To ensure that
Wij is in a biophysically reasonable range, we
need a relation between it and the sizes of both excitatory and
inhibitory PSPs. To derive such a relation, we use Eq. 5a to
compute, as a function of Wij, the peak voltage in response to a single presynaptic action potential. For definiteness we consider a neuron whose resting membrane potential is the nominal one, Vr (which implies that
Îa = 0), and assume that the neuron has
not produced an action potential for a sufficiently long time that both
K and
KCa have decayed to zero.
Under these conditions, Eq. 5a for the membrane potential can be combined with Eq. 6 for the synaptic current to yield
Letting Vi Vi Vr, linearizing the
resulting equation, and assuming the synaptic drive is small, we find
that Vi evolves according to
If neuron j fires at time t = 0, the
fraction of open channels, sij, jumps
instantaneously from 0 to rs and subsequently decays exponentially: sij(t < 0) = 0, sij(t 0) = rs
exp(t/s) (recall that
s is the decay time of the open channels and
rs is the use parameters; see Eq. 7).
Combining this time dependence for sij with the
initial condition Vi (t = 0) = 0, the above equation has the solution
This expression implies that
|Vi(t)| is maximum when
t = ln
(cell/s)/(s1 cell1) t0.
Defining VPSPij Vi(t0), it is
straightforward to show that
|
(14)
|
In our simulations, rather than specifying the size of
Wij directly, we specify the sizes of the
excitatory and inhibitory PSPs
(VPSPij with the type of neuron
j chosen appropriately) and then use Eq. 14 to
determine Wij.
PARAMETERS.
Table 1 contains a list of all parameters
used in our simulations. Parameters followed by an asterisk are the
ones we vary; for those, a range or set of values is given. Parameters
not followed by an asterisk remain constant throughout the simulations.
The grouping into single-cell, synaptic, and network parameters is in
most cases clear; the exceptions are VEPSP and
VIPSP, the nominal amplitudes of the excitatory
and inhibitory postsynaptic potentials. These are usually thought of as
synaptic or single-cell properties, but in our simulations we use them
only to determine the connectivity matrix, Wij,
via Eq. 14. Specifically,
VPSPij = VEPSP when neuron j is excitatory and
VPSPij = VIPSP when neuron j is inhibitory.
Thus they are listed under network parameters.
We performed simulations with two networks, which we denote
network A and network B. The former was chosen to
simulate networks of cortical neurons, the latter to simulate networks
of cultured mouse spinal cord neurons in which we did experiments
relevant to the theoretical predictions made here (Latham et al.
2000). They differ primarily in their connectivity and
postsynaptic potential amplitude, network A having high,
infinite range connectivity and small PSPs and network B
having low, local connectivity and large PSPs. The parameters specific
to each network are given in Table 2,
which contain all the parameters marked with an asterisk from Table 1.
A dash in Table 2 indicates parameters that are varied during the
course of the simulations.
Most of the parameters in Tables 1 and 2 are either dimensionless or
have physical units. There are, however, three normalized parameters:
Îmax, which has units of voltage, and
K and
KCa, which are both
dimensionless. To convert these to their correct physical units (Amps
for current and Siemans for conductance) divide by the cell membrane resistance.
Simulations were performed with a fourth-order Runge-Kutta integration
scheme. In all simulations reported we used a time step of 1 ms. We
occasionally checked that this was short enough by decreasing the time
step to 0.5 ms, and in no cases did we see a change in our results.
Choice of simulation parameters
In choosing parameters we attempted to stay close to values
observed in two different systems: networks of cortical neurons, and
networks of cultured mouse spinal cord neurons in which we did
experiments relevant to the theoretical predictions made here (Latham et al. 2000). These two systems differ primarily
in their connectivity and postsynaptic potential amplitude, the former having high connectivity and small PSPs, the latter low connectivity and large PSPs. Thus most parameters were the same for the two model
networks used in our simulations, and these were relatively standard;
the parameters that differed pertained to connectivity and PSP
amplitude, mirroring the differences between cortical and cultured
spinal cord networks.
As is well-known, connectivity in cortex is high; each neuron connects
to approximately 7,000 others (Braitenberg and Schüz 1991). The size of a local circuit, assuming such a thing
exists, is more difficult to determine. However, given that axonal
spread is measured in hundreds of micrometers and the density of
neurons is ~105/mm3 (Braitenberg and
Schüz 1991), a local circuit on the order of 100,000 neurons is a reasonable estimate. Unfortunately, we do not yet have the
computational power to model such large networks. Consequently, we
considered networks of only 10,000 neurons, each of which made ~1,000
connections. To make up for the reduced connectivity in our model
network we increased the size and frequency of the postsynaptic
potentials: instead of using the more realistic numbers of 0.2 mV for
EPSPs (Komatsu et al. 1988; Matsumura et al.
1996) and 1 mV for IPSPs (Tamás et al.
1998), we typically used 1.0 and 1.5 mV. The frequency of
PSPs was increased in our model by ignoring failures, so
neurotransmitter was released every time an action potential occurred.
Unlike cortex and spinal cord in vivo, cultured spinal cord networks
are intrinsically localized. However, connectivity in our particular
preparation has not been characterized. We thus made estimates as
follows. The number of connections a neuron makes is
APA, where A is the
area of the axonal arborization, is the number of neurons per area,
and PA is the probability that a neuron connects
to another neuron within its axonal arborization. We measured axonal
arborization from neurons visualized after intracellular injection with
horseradish peroxidase (Neale et al. 1978), and found it
to be 2.0 ± 0.86 mm2 (mean ± SD,
n = 12 neurons). The density of our cultures, , was
23.7 ± 10.2 neurons/mm2 (n = 10 culture dishes). Thus a single neuron is capable of connecting to
~475 others. Finally, in paired recordings (Nelson et al.
1981), approximately one-half the cell pairs tested were
connected. This gives a connectivity of ~240 connections/neuron. We
used 200 connections/neuron in most of our simulations, although we
explored a range of values.
The size of the postsynaptic potentials is large in cultured spinal
cord neurons; on the order of 3-10 mV for EPSPs and approximately 6
mV for IPSPs (Nelson et al. 1981, 1983). In contrast to cortex, transmission failures
are not observed; this is because a presynaptic neuron makes a large
number of connections (~100) on each of its postsynaptic targets
(Nelson et al. 1983; Pun et al. 1986). In our simulations we typically used EPSP and IPSP amplitudes of 4 and 6
mV, respectively, and, for simplicity, we did not include variation.
The parameters given in Tables 1 and 2 were the ones we used in the
bulk of our simulations. In APPENDIX A we studied the
effects of varying these parameters. In no cases did changing parameters produce network behavior that was inconsistent with our
model. This suggests that our model made robust predictions and that
our simulation results were not due simply to a particular choice of parameters.
Reduced network model
The reduced network model we use in our analysis is based on the
Wilson and Cowan equations (Wilson and Cowan 1972), in
which the dynamical variables are the mean excitatory and inhibitory firing rates. Those equations are given in generic form in Eq. 16 of RESULTS. Here we write down a particular
realization of Wilson and Cowan-type equations that we augment by
including spike-frequency adaptation.
Although it would be desirable to start with the equations describing
the full network simulations and derive from those a set of mean firing
rate equations, this is essentially impossible to do analytically
except in highly idealized cases (van Vreeswijk and Sompolinsky
1998), and extremely difficult numerically. Thus we propose a
set of reduced equations designed to capture qualitative, but not
quantitative, features of large-scale networks. Denoting 
E and I
as the mean excitatory and inhibitory firing rates (averaged over
neurons), and GE and GI
as the mean level of spike-frequency adaptation of the excitatory and
inhibitory populations, we let these variables evolve according to
![&tgr;<SUB><IT>E</IT></SUB> <FR><NU><IT>d<A><AC>&ngr;</AC><AC>&cjs1171;</AC></A></IT><SUB><IT>E</IT></SUB></NU><DE><IT>d</IT><IT>t</IT></DE></FR><IT>=</IT><IT>g</IT>[<IT>A</IT>(<IT>C<SUB>EE</SUB></IT><IT><A><AC>&ngr;</AC><AC>&cjs1171;</AC></A></IT><SUB><IT>E</IT></SUB><IT>−</IT><IT>C<SUB>EI</SUB></IT><IT><A><AC>&ngr;</AC><AC>&cjs1171;</AC></A></IT><SUB><IT>I</IT></SUB><IT>+&thgr;<SUB>max</SUB>−</IT><IT>G<SUB>E</SUB></IT>)<IT>/</IT>(<IT>1+<A><AC>&ngr;</AC><AC>&cjs1171;</AC></A></IT><SUB><IT>I</IT></SUB>)]<IT>−<A><AC>&ngr;</AC><AC>&cjs1171;</AC></A></IT><SUB><IT>E</IT></SUB>](/content/vol83/issue2/fulltext/808/img035.gif)
|
(15a)
|
![&tgr;<SUB><IT>I</IT></SUB> <FR><NU><IT>d<A><AC>&ngr;</AC><AC>&cjs1171;</AC></A></IT><SUB><IT>I</IT></SUB></NU><DE><IT>d</IT><IT>t</IT></DE></FR><IT>=</IT><IT>g</IT>[<IT>A</IT>(<IT>C<SUB>IE</SUB></IT><IT><A><AC>&ngr;</AC><AC>&cjs1171;</AC></A></IT><SUB><IT>E</IT></SUB><IT>−</IT><IT>C<SUB>II</SUB></IT><IT><A><AC>&ngr;</AC><AC>&cjs1171;</AC></A></IT><SUB><IT>I</IT></SUB><IT>+&thgr;<SUB>max</SUB>−</IT><IT>G<SUB>I</SUB></IT>)<IT>/</IT>(<IT>1+<A><AC>&ngr;</AC><AC>&cjs1171;</AC></A></IT><SUB><IT>I</IT></SUB>)]<IT>−<A><AC>&ngr;</AC><AC>&cjs1171;</AC></A></IT><SUB><IT>I</IT></SUB>](/content/vol83/issue2/fulltext/808/img036.gif)
|
(15b)
|
|
(15c)
|
|
(15d)
|
where E and
I are the time scales for relaxation to a
firing rate equilibrium, SFA is the
characteristic time scale for changes in the level of spike-frequency
adaptation, A is the overall amplitude of the gain
functions, the CLM, L, M = E, I, correspond to
connectivity among the excitatory (E) and inhibitory
(I) populations, G corresponds to the coupling between firing rate and spike-frequency adaptation, max
controls network excitability, and g(x) is a
thresholding function that saturates at 
max, the maximum
firing rate of the neurons,
The mean levels of spike-frequency adaptation,
GE and GI, are related
loosely to gKCa in the network simulations, G is related to
KCa, and
max is related to Îmax.
The qualitative features that these equations capture are 1)
the gain functions, g(...), are generally increasing
functions of excitatory firing rate and decreasing functions of
inhibitory firing rate [the CLM are all
positive and g(x) is a nondecreasing function of
x], 2) increasing the firing rate increases
spike-frequency adaptation (G is positive), 3)
increasing the level of spike-frequency adaptation reduces firing
rates, 4) the gain functions are approximately threshold-linear, consistent with the numerical gain functions in Fig.
1A, and 5) the gain is reduced by a divisive
term proportional to the inhibitory firing rate. The divisive term,
which has been proposed as a mechanism for gain control in cortical
neurons (Carandini and Heeger 1994; Carandini et
al. 1997; Heeger 1992; Nelson
1994), was included for two reasons. First, it was observed in
our simulations: Fig. 1A in APPENDIX B shows a
pronounced drop in gain as the inhibitory firing rate increases.
Second, curved nullclines are essential for the existence of the
saddle-node bifurcation that produces bursting (Fig. 5), and the
divisive term increases nullcline curvature, thus making bursting more robust.
We are interested in understanding the dynamics represented by
Eq. 15 in the regime in which E
and I are on the order of the membrane time
constant, ~10 ms, and SFA is on the order
of the spike-frequency adaptation time, which can be as high as several
seconds. Thus E,
I SFA, and
we can use the fast/slow dissection proposed by Rinzel
(1987) to analyze the dynamics. With this approach, the full
four-dimensional system is reduced to two subsystems: a fast one
corresponding to the population firing rates,
E and
I, and a slow one corresponding to the
level of spike-frequency adaptation, GE and
GI. Although this represents a major
simplification, it still leaves us with three main behaviors: steady
firing, oscillations, and bursting. The first two can be readily
understood in terms of the two-dimensional firing rate dynamics with
GE and GI held fixed;
only the third, bursting, requires the interaction of all four
variables. We thus briefly discuss a method for its analysis.
The idea behind the fast/slow dissection is that
E and I
relax rapidly to an attractor (either a steady state or oscillations)
compared with the time scale over which GE and
GI change. Motion in the
GE GI plane is
then determined be Eqs. 15c and 15d with
E and I
evaluated on their attractor. In the purely bursting regime, for
attainable levels of spike-frequency adaptation the network is either
silent or fires steadily; the attractor is a fixed point for given
values of GE and GI. In
this steady-state regime, the mean inhibitory firing rate,
I, is a single-valued function of the mean excitatory rate, E (see Fig.
3B). The dynamics associated with Eq. 15 can thus
be reduced to three dimensions, and the important features of the
dynamics are captured by a three-dimensional bifurcation diagram that
consists of a sheet whose height above the GE GI plane represents the value of
E. Bifurcations (sudden changes in
firing rate) occur at folds in the sheet.
For a particular set of parameters, the trajectory in the
GE GI plane is
typically a closed curve that exhibits sudden changes of direction at
bifurcation points. If the trajectory collapses onto, or almost onto, a
single curve, we may make a further, approximate, reduction to a
two-dimensional bifurcation diagram of
E versus GE. For
the parameters considered in RESULTS (Fig. 6), the curve in
the GE GI plane
collapsed almost to a line, resulting in the two-dimensional
bifurcation diagrams shown in Fig. 6, C-E. For each of
these diagrams the lines were found numerically using least-squares
regression. The equations for the lines were as follows: Fig.
6C (max = 0), GI = 0.00 + 0.64GE,
R2 = 1.0000; Fig. 6D
(max = 0.25), GI = 0.01 + 0.65GE, R2 = 0.9996; Fig. 6E (max = 0.52),
GI = 0.03 + 0.65GE,
R2 = 0.9999.
| |
RESULTS |
In this section we investigate how single neuron and network
properties affect intrinsic firing patterns. We use as our model system
large, isolated networks with random, sparse connectivity and develop a
theory that describes the firing patterns in such networks in the low
firing rate regime. We then test specific predictions of this theory by
performing simulations with large networks of spiking neurons.
Robustness of the theory is verified in APPENDIX A, where
we examine a broad range of single neuron and network parameters.
Theory
The theoretical development is divided into three parts:
1) analysis of networks that do not exhibit adaptation,
2) analysis of the more realistic, and more biologically
relevant, case of networks that do exhibit adaptation, and
3) analysis of the stability of low firing rate equilibria.
LOW FIRING RATE EQUILIBRIA IN SPARSE, RANDOMLY CONNECTED
NETWORKS WITHOUT ADAPTATION.
To describe the dynamics of large, sparse, randomly connected networks,
we start with the model proposed by Wilson and Cowan (1972). We use as our dynamical variables the mean excitatory <
TOP
ABSTRACT
INTRODUCTION
![<IT>Z</IT><SUB><IT>T<SUB>j</SUB></IT></SUB><IT>=</IT><LIM><OP>∫</OP><LL><IT>0</IT></LL><UL><IT>2&pgr;</IT></UL></LIM> <FR><NU><IT>d&thgr;</IT></NU><DE><IT>&pgr;</IT></DE></FR> <LIM><OP>∫</OP><LL><IT>0</IT></LL><UL><IT>1</IT></UL></LIM> <IT>rdr</IT><IT> exp</IT>(−<IT>r</IT><SUP><IT>2</IT></SUP><IT>/2&sfgr;</IT><SUP><IT>2</IT></SUP><SUB><IT>T<SUB>j</SUB></IT></SUB>)<IT>=2&sfgr;</IT><SUP><IT>2</IT></SUP><SUB><IT>T<SUB>j</SUB></IT></SUB>[<IT>1−exp</IT>(−<IT>1/2&sfgr;</IT><SUP><IT>2</IT></SUP><SUB><IT>T<SUB>j</SUB></IT></SUB>)]](/content/vol83/issue2/fulltext/808/img029.gif)
