INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B REFERENCES

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

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B REFERENCES

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, N_E of which are excitatory and N_I of which are inhibitory. The time evolution equation for the membrane potential of neuron i, V_i, may be written

(1)

where C_cell is the cell capacitance, I_spike consists of the currents responsible for generating action potentials, I_fAHP and I_sAHP are fast and slow afterhyperpolarization currents, respectively, I_syn represents the synaptic current, and I_a 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.

(2)

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.

(3a)

(3b)

KCa,i

(4a)

(4b)

K,

KCa,

KCa

KCa

KCa

KCa,i

K,

KCa,

KCa

KCa

KCa,i

K,

KCa,

(5a)

(5b)

(5c)

_K

KCa

_K

KCa

(6)

(7)

(8)

(9a)

(9b)

_i

,i

(10a)

(10b)

_i

,i

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, W_ij. 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 P_ij. 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, W_ij is set to zero.

(11a)

(11b)

(12a)

(12b)

(13)

_K

KCa

1

1

(14)

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 V_EPSP and V_IPSP, 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, W_ij, via Eq. 14. Specifically, V_PSPij = V_EPSP when neuron j is excitatory and V_PSPij = V_IPSP when neuron j is inhibitory. Thus they are listed under network parameters.

_K

KCa

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 ~10⁵/mm³ (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 AP_A, where A is the area of the axonal arborization,  is the number of neurons per area, and P_A 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 mm² (mean ± SD, n = 12 neurons). The density of our cultures, , was 23.7 ± 10.2 neurons/mm² (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 G_E and G_I as the mean level of spike-frequency adaptation of the excitatory and inhibitory populations, we let these variables evolve according to

(15a)

(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 C_LM, 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, G_E and G_I, are related loosely to g_KCa 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 C_LM 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, G_E and G_I. 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 G_E and G_I 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 G_E and G_I change. Motion in the G_E  G_I 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 G_E and G_I. 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 G_E  G_I 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 G_E  G_I 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 G_E. For the parameters considered in RESULTS (Fig. 6), the curve in the G_E  G_I 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), G_I = 0.00 + 0.64G_E, R² = 1.0000; Fig. 6D (_max = 0.25), G_I = 0.01 + 0.65G_E, R² = 0.9996; Fig. 6E (_max = 0.52), G_I = 0.03 + 0.65G_E, R² = 0.9999.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B REFERENCES

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 and inhibitory firing rates, denoted _E and _I, respectively. In terms of these variables, the Wilson and Cowan model can be cast in the form

(16a)

(16b)

where _E and _I are time constants that determine how fast the network relaxes to its equilibria and _E and _I are the excitatory and inhibitory gain functions. The gain functions determine the network firing rates at one instant of time given the firing rates at an earlier time. Specifically, if the average excitatory and inhibitory firing rates at time t are _E and _I, then _E and _I are the average excitatory and inhibitory firing rates at times t + _E and t + _I, respectively. In the original Wilson and Cowan equations the gain functions consisted of two terms, one associated with the neurons' refractory period and one with the "subpopulation response function" (Wilson and Cowan 1972). Because we are interested in low firing rates, we have ignored the term associated with the refractory period.

(17a)

(17b)

View larger version (20K): [in this window] [in a new window]   Fig. 2. Schematic of the excitatory and inhibitory gain functions in the 3 regimes discussed in the main text. A: excitatory gain functions. In regime 1, E(E, 0) is zero when E = 0 (neurons do not fire without input), and its slope at E = 0 is also zero (a single action potential cannot make a neuron fire, and thus cannot produce any activity in a silent network). In regime 2 the neurons still cannot fire without input, so E(0, 0) is again zero. However, because a single action potential can cause cells to fire and thus ignite network activity, the slope at E = 0 is greater than one. The condition for a slope greater than one is that at least 1/KEE of the excitatory neurons have their resting membrane potential within one excitatory postsynaptic potential (EPSP) of threshold, where KEE is the mean number of excitatory connections made by an excitatory neuron. Because KEE is large, this could easily occur. In regime 3 some of the neurons fire without input, so E(0, 0) > 0. B: Inhibitory gain functions. In regimes 1 and 2, I(0, I) = 0, whereas in regime 3, I(0, 0) > 0.

View larger version (25K): [in this window] [in a new window]   Fig. 3. Gain functions and nullclines when none of the cells have their resting membrane potential within one EPSP of threshold (regime 1): A: excitatory gain curves, E(E, I), vs. E for various values of I. Only the curve with I = 0 is labeled; the other curves correspond to increasingly larger values of I. Each intersection between a gain curve and the line E = E corresponds to a point on the excitatory nullcline, as indicated by the connecting dashed lines. B: inhibitory gain curve, I(E, I), vs. I for various values of E. Again, only the curve with E = 0 is labeled (the straight line at I = 0); the other curves corresponds to increasingly larger values of E. Again, each intersection between a gain curve and the line I = I corresponds to a point on the excitatory nullcline, as indicated by the connecting dashed lines. C: the associated excitatory (gray) and inhibitory (black) nullclines. These nullclines separate regimes where the time derivatives of the firing rates are positive from regimes where they are negative: dE/dt is positive below the excitatory nullcline and negative above it, whereas dI/dt is positive to the right of the inhibitory nullcline and negative to its left. This rule is illustrated by the horizontal and vertical arrows, which indicate the signs of dE/dt and dI/dt, respectively.

View larger version (26K): [in this window] [in a new window]   Fig. 4. Nullclines in the 3 regimes described in the main text. Excitatory and inhibitory nullclines are shown with gray and black lines, respectively. As in Fig. 3, dE/dt is positive below the excitatory nullcline and negative above it, whereas dI/dt is positive to the right of the inhibitory nullcline and negative to its left. S, M, and U indicate stable, metastable, and unstable equilibria, respectively. A: regime 1, none of the cells have their resting membrane potential within one EPSP of threshold. There is a stable equilibrium at zero firing rate and a metastable one at high firing rate, but the low firing rate equilibrium is unstable. B: regime 1, with strong curvature introduced in the inhibitory nullcline to create a metastable, low firing rate state. C: regime 1, in the very high connectivity limit where the nullclines are straight. In this limit 3 properties conspire to eliminate the possibility of a low firing rate equilibrium: the inhibitory nullcline is tied to the origin, there is a gap between the origin and the unstable branch of the excitatory nullcline, and the nullclines are straight. D: regime 2, some cells are within one EPSP of threshold. A metastable state can exist. This state persists even in the high connectivity limit, where the nullclines become straight. E: regime 3, endogenously active cells are present. A single, globally attracting equilibrium can exist. F: same as E except in the very high connectivity regime.

EFFECT OF ADAPTATION ON FIRING PATTERNS. So far we have assumed that the gain curves, and thus the nullclines, are static. In fact, this is not the case for real neurons. Because of spike-frequency and/or synaptic adaptation, nullclines are dynamic objects that depend on the history of activity, not just the activity at a single point in time. To incorporate this history into our analysis, we assume that neurons adapt slowly compared with the time it takes a network to reach equilibrium. Consequently, at any instant in time the equilibrium firing rates occur at the intersection of the excitatory and inhibitory nullclines, just as they did above. However, the equilibria are not fixed: adaptation causes slow changes in the nullclines, and thus in the equilibrium firing rates.

View larger version (16K): [in this window] [in a new window]   Fig. 5. Excitatory and inhibitory nullclines at fixed levels of adaptation in a regime where adaptation causes bursting. For clarity, only the parts of the nullclines at low firing rate are shown. Because the level of adaptation changes slowly, we refer to the intersections of the excitatory and inhibitory nullclines as equilibria. These "equilibria," which in fact change with time, are marked by black dots. A: the level of adaptation is minimum, and there is only 1 firing rate equilibrium; the equilibrium near zero has just disappeared via a saddle-node bifurcation. B: the level of adaptation has increased, leading to a bistable state. Barring large fluctuations, the network will stay near the higher firing rate state. C: the level of adaptation is large enough that the equilibrium at higher firing rate is eliminated, again via a saddle-node bifurcation. The resulting crash to zero causes the level of adaptation to begin to decrease. With time, this results in a shift to the low firing rate equilibrium in B. With a further decrease in adaptation, the low firing rate equilibrium disappears (A), and the firing rate jumps to a higher value. At this point the cycle repeats.

View larger version (22K): [in this window] [in a new window]   Fig. 6. Particular realization of the Wilson and Cowan equations, Eq. 15 of METHODS, showing transitions from steady firing to bursting to silence. A: equilibrium firing rate of the excitatory population as a function of max, which regulates network excitability: max > 0 corresponds to the existence of endogenously active cells, max < 0 corresponds to their absence (see text and METHODS). When max > 0.51 there is a steady, low firing rate equilibrium, indicated by the single solid line to the right of the gray region. As max decreases, there is a transition to bursting. In the bursting regime (0 < max < 0.51) the firing rate periodically jumps between 0 Hz and the solid black line (the mean firing rate averaged over a burst) above the gray region. When max < 0, corresponding to the absence of endogenously active cells, the network crashes to zero firing rate and stays there. B: mean excitatory firing rate (E, black line) and mean level of excitatory spike-frequency adaptation (GE, gray line) vs. time in the bursting regime, max = 0.25. Note that the level of spike-frequency adaption, GE, increases when the network is active and decreases when it is silent. C-E: approximate 2-dimensional bifurcation diagrams (see METHODS). Black solid curves in each plot denote the stable branches of the excitatory firing rate equilibrium; black dashed curve that connects them denotes the unstable one. Thick gray line is the spike-frequency adaptation nullcline associated with Eq. 15c. Thin dashed gray line is the trajectory. C: max = 0, so there is a steady-state equilibrium at zero firing rate. A trajectory is shown that starts near E = 6, GE = 0, traverses the bifurcation diagram, and ends at the point marked "S" (for stable) at E = GE = 0. D: max = 0.25, placing the network in the bursting regime. The trajectory in this regime, which cycles indefinitely, reveals a relaxation oscillator; we refer to the resulting dynamics as bursting because individual neurons fire repetitively during the active phase. E: max = 0.52, placing the network just inside the steady firing rate regime. A trajectory is shown that starts near E = 0, GE > 0, traverses the bifurcation diagram, and ends at the point marked "S." The transition between bursting and steady firing occurs at the point predicted by the bifurcation diagram: the spike-frequency adaptation nullcline (the heavy gray line) passes onto the stable branch of the excitatory firing rate equilibrium at max  0.52. Parameters (see Eq. 15): E = I = 10 ms, SFA = 2,000 ms, A = 20, CEE = 0.9, CEI = 1.0, CIE = 1.0, CII = 1.4, G = 0.2, and max =  (chosen large because the firing rate is far from saturation). Plots were made using XPP, developed by G. B. Ermentrout.

STABILITY OF LOW FIRING RATE EQUILIBRIA. Although endogenously active cells are necessary for the existence of a stable, low firing rate equilibrium, such cells do not guarantee either of these properties. As parameters change, a low firing rate equilibrium can become unstable via a Hopf bifurcation (Marsden and McCracken 1976), or it can be driven to high firing rate. In this section we investigate stability, then discuss conditions under which a network is both stable and fires at low rate.

(18)

(19)

SUMMARY OF THEORETICAL ANALYSIS. The tool we used to analyze network behavior was a simplified firing rate model, in which the dynamics of networks containing greater than 10⁴ neurons was reduced to a small number of equations describing average quantities. Rather than using a particular firing rate model, we used geometrical arguments to derive generic network behavior. Our primary assumptions were that, on average, excitatory input to a neuron increases its firing rate, inhibitory input decreases it, and neurons exhibit spike-frequency adaptation. With these assumptions, we were able to show that 1) endogenous activity is necessary for low firing rates, 2) decreasing the fraction of endogenously active cells causes a network to burst, and 3) firing patterns, especially transitions between steady firing and oscillations, are dependent on network connectivity. A specific realization of a simplified firing rate model corroborated these predictions.

Simulations

To test the predictions made in the previous section, we performed simulations with large networks of spiking model neurons, as described in METHODS. We examined 1) the effect of the distribution of applied depolarizing currents, I_a, on firing patterns and 2) how relative connection strengths among excitatory and inhibitory populations affect network behavior.

We used two networks, denoted A and B. These networks differed primarily in their connectivity and the sizes of the PSPs of their constituent neurons. Network A was designed to simulate cortical networks. It contained 8,000 excitatory and 2,000 inhibitory neurons, and each neuron connected, on average, to 1,000 others. Connectivity was infinite range (meaning the probability of 2 neurons connecting did not depend on the distance between them; see METHODS), and PSPs were on the order of 1 mV. Network B was chosen to be consistent with our experiments (Latham et al. 2000) and to test whether the predictions made for infinite range connectivity apply also to local connectivity, for which the connection probability depends on distance between neurons. It had fewer connections per neuron than Network A (200 instead of 1,000), local rather than infinite range connectivity, larger PSPs (on the order of 5 mV rather than 1 mV), and a higher fraction of inhibitory neurons (30% instead of 20%). A complete list of the parameters for each network is given in Tables 1 and 2 of METHODS.

Networks A and B behaved similarly, so in the bulk of this section we concentrate on Network A. At the end of the section we summarize the results of Network B.

ENDOGENOUS ACTIVITY. For networks without spike-frequency adaptation, the theoretical analysis presented above led to the following predictions. If endogenously active cells are present, then a stable, low firing rate equilibrium is possible. If endogenously active cells are absent, a low firing rate equilibrium is still possible. However, it is metastable; the network will eventually decay to a zero firing rate state. In this case there is a further distinction between networks in which some cells have their resting membrane potential within one EPSP of threshold and networks in which none do. This distinction has primarily qualitative bearing on firing patterns. In networks containing cells that are one EPSP from threshold, a single action potential can ignite activity in a silent network, so finite firing rate states tend to be long lived. In networks with no cells within one EPSP of threshold, finite firing rate states tend to have shorter lifetimes.

View larger version (18K): [in this window] [in a new window]   Fig. 7. Mean excitatory firing rate () vs. Îmax for Network A. (The inhibitory rate, which is not shown, is ~2.5 times higher for all data points.) The gray region in B indicates bursting networks; for such networks, the top and bottom squares show the firing rates in the active and silent phases of the burst, respectively. Squares at zero that are not in the gray region indicate networks that crashed to zero and remained there. Each firing pattern has associated with it a set of nullclines with a particular underlying structure (see Fig. 4). Firing patterns are matched to generic examples of their associated nullclines with arrows. Bursting networks alternate between a set of nullclines with a stable, low firing rate equilibrium, and a set of nullclines in which the only equilibrium is at zero firing rate. As discussed in the main text (see especially Fig. 5B), between these 2 extremes is a bistable state supporting both steady firing and silence. A: no spike-frequency adaptation. The networks either fire steadily for 100 s, the length of the simulations, or crash to zero firing rate. Some equilibria with Îmax < Îthresh (no endogenously active cells) are metastable but relatively long-lived. B: the cells exhibit spike-frequency adaptation. This abolishes the long-lived metastable states and allows bursting. For both A and B the parameters were taken from Network A of Table 2 with BE = 1.4 and BI = 1.1; the difference between the 2 networks lies in the amount, KCa, the normalized potassium conductance increased on each spike. In A, KCa = 0; in B, KCa = 0.01. For these network parameters, ÎE1  3.715 (determined numerically) and Îthresh = 3.75 (determined analytically; see the section DISTRIBUTION OF APPLIED CURRENTS in METHODS). To convert Îmax, ÎE1, and Îthresh from their normalized units, mV, to nA, divide by the cell membrane resistance in M.

View larger version (34K): [in this window] [in a new window]   Fig. 8. Firing patterns in the bursting regime. Top: raster plots of 200 randomly chosen neurons in the burst regime. The 1st 40 (1-40) are inhibitory and the last 160 are excitatory. Bottom: voltage trace for neuron 40. During the burst, voltage fluctuations are PSP driven; between bursts, the whole network is silent so the voltage decays to the resting membrane potential, 50 mV. Because there is no noise in these simulations, between bursts the voltage is constant. Parameters are taken from Fig. 7 with Îmax = 3.82.

CONNECTIVITY. The theoretical predictions concerning the effect of coupling on stability and firing rate were 1) decreasing I-I and increasing E-E coupling both lead to oscillations in firing rate and 2) increasing same-type coupling (E-E or I-I) causes mean firing rates to go up, whereas increasing opposite-type coupling (E-I or I-E) causes mean firing rates to go down.

View larger version (28K): [in this window] [in a new window]   Fig. 9. Mean firing rate vs. time for excitatory () and inhibitory (- - -) populations at different value of the bias parameters, BE and BI. A and D: because the I-I coupling is low (small BI), these networks oscillate. B, C, and E: with increasing I-I coupling (increasing BI), oscillations are eliminated and average firing rates rise. The small fluctuations in these plots are due to the finite size of the network. F: larger I-I coupling further raises average firing rates and produces bursting. The gray line on this plot is the normalized slow potassium conductance averaged over excitatory neurons (see Eq. 5c). This quantity corresponds to the mean level of spike-frequency adaptation, so it increases when the neurons are firing and decreases when they are silent. (To clearly show the bursting, the time scale in this panel is different from the other ones and the inhibitory firing rate, which was ~30% lower, is not shown.) The difference between oscillations (A and D) and bursting (F) is that in the former the rise and fall of the firing rate is slow compared with the period, whereas in the latter it is fast. Parameters for all panels were taken from Network A of Table 2 with Îmax = 5.0 and KCa = 0.01. Firing rates were computed by convolving the spike trains with a Gaussian of width 10 ms.

KCa,E

KCa,E

1

iE

KCa,i

KCa,E

B,

KCa,E

KCa

LOCAL CONNECTIVITY. The above simulations were repeated with parameters corresponding to Network B, in which the connectivity is local rather than infinite range (meaning neurons that are close together are more likely to connect than those that are far apart), the number of connections per neuron is smaller than in Network A, and the PSPs are larger (see Tables 1 and 2). The results are summarized in Figs. 10 and 11.

View larger version (17K): [in this window] [in a new window]   Fig. 10. Mean excitatory firing rate () vs. Îmax for Network B. (The inhibitory rate, which is not shown, is a factor of ~2.5 higher for all data points.) See Fig. 7 caption for details. For both A and B the parameters were taken from Network B of Table 2 with BE = 1.4 and BI = 1.0; the difference between the 2 networks is that in A, KCa = 0, whereas in B, KCa = 0.08. For these network parameters, ÎE1  3.230 (determined numerically) and Îthresh = 3.75 (determined analytically).

View larger version (30K): [in this window] [in a new window]   Fig. 11. Mean firing rate vs. time for excitatory () and inhibitory (- - -) populations at different value of the bias parameters, BE and BI. In contrast to Network A, Fig. 9, this network (Network B) did not oscillate at any of the parameters shown, although oscillations were observed at lower values of the I-I coupling (smaller BI). A, D, and E: steady firing. B, C, and F: bursting. As in Fig. 9F, the gray line is the mean value of the normalized slow potassium conductance averaged over excitatory neurons. Note the different time scales than in A, D, and E, and that the inhibitory firing rate is not shown. The parameters for all panels were taken from Network B of Table 2 (local connectivity) with Îmax = 5.0 and KCa = 0.08. The firing rates were computed by convolving the spike trains with a Gaussian of width 10 ms.

KCa,E

SUMMARY OF SIMULATION RESULTS. We performed simulations with two very different networks: an infinite range, high connectivity network with small PSPs, and a local, low connectivity network with large PSPs. We also examined a broad range of parameters, including the distribution of applied currents, spike-frequency adaptation, and network connectivity. The results of all simulations were as predicted: for networks without spike-frequency adaptation, we observed transitions from silence to steady firing as the fraction of endogenously active cell increased. For networks with spike-frequency adaptation, there was an intermediate regime in which the networks burst. Finally, when the I-I coupling was strengthened, we observed transitions from oscillations to steady firing to bursting and an increase in firing rate.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX A APPENDIX B REFERENCES

We have investigated, both theoretically and through simulations with spiking model neurons, the intrinsic dynamics in large networks of neurons. The goal of this work was twofold: 1) to understand how dynamic interactions between excitatory and inhibitory neurons lead to stable, low firing rates, such as those widely observed in the mammalian central nervous system, and 2) to determine the conditions under which networks switch from steady firing to bursting.

An understanding of the mechanism for generating stable, low firing rates in large, isolated neuronal networks has been elusive (Abeles 1991; Amit and Treves 1989; Buhmann 1989; Treves and Amit 1989). The elusiveness stems from the difficulty in controlling the powerful recurrent excitation that exists in such networks. To compensate for this recurrent excitation, inhibitory feedback is required. To stabilize low firing rates, that feedback must be strong enough that the inhibitory firing rate is more sensitive to input from excitatory cells than the excitatory firing rate. Or, in the language of dynamics, the inhibitory nullcline in average firing rate space must be steeper than the excitatory one at the equilibrium (e.g., Fig. 4E). We found that for isolated networks, the above condition occurs robustly at low firing rate only when endogenously active cells are present; without such cells networks are either silent or fire at high rate. This led to the following prediction: if low firing rates are observed in an isolated network, then that network must contain endogenously active cells. This prediction was confirmed by large-scale network simulations, which included the exploration of a broad range of parameters to ensure that the simulations were robust (APPENDIX A). It was also corroborated by experiments in cultured neuronal networks, as described in the accompanying paper (Latham et al. 2000).

Although endogenously active cells are necessary for the existence of a low firing rate equilibrium, they do not guarantee its stability. In particular, we found that a high level of spike-frequency adaptation leads to bursting: repetitive firing can introduce a hyperpolarizing current sufficient to temporarily eliminate endogenously active cells, resulting in a crash to zero firing rate; after the cells stop firing, the hyperpolarizing current decays and firing resumes. (Interestingly, our analysis implies that synaptic adaptation, because it does not affect the fraction of endogenously active cells, cannot produce bursting.) The probability of bursting is highest if there are few endogenously active cells to begin with, and for that reason the fraction of such cells plays a key role in determining firing patterns. Specifically, networks with no endogenously active cells are typically silent; at some finite fraction of endogenously active cells there is a transition to bursting; and at an even higher fraction there is a second transition to steady firing at low rate. This scenario was also confirmed by network simulations, and it was corroborated by experiments in neuronal cell culture (Latham et al. 2000).

Although the fraction of endogenously active cells plays the dominant role in determining the primary firing patterns (silence, bursting, or steady firing), another parameter, network connectivity, influences secondary features of the firing patterns. Specifically, when inhibitory-inhibitory coupling becomes too weak or excitatory-excitatory coupling becomes too weak or excitatory-excitatory coupling becomes too strong, a low firing rate equilibrium can be destabilized via a Hopf bifurcation. This leads to oscillations in firing rate that can occur in both the "steady" firing and the bursting regimes. In the latter case, the oscillations occur during the active phase of the burst.

These theoretical findings were based on the analysis of isolated networks with random, infinite-range connectivity. However, theoretical considerations indicate that they apply to networks that receive external input, and simulations indicate that they are valid for networks with local connectivity, in which the connection probability is a decreasing function of distance between neurons. The validity of the model for local connectivity suggests that these findings will hold for the more structured architecture that exists in the cortex.

Implications

Two firing patterns that are ubiquitous in the mammalian CNS are steady firing at low rates and rhythmic bursting. We have shown that, in isolated networks, both firing patterns are controlled largely by a single parameter, the fraction of endogenously active cells. As long as the fraction of such cells is above a threshold, steady firing is possible; below that threshold the network bursts. The threshold depends on the degree of spike-frequency adaptation, a property of neurons that has been shown to be modulatable [e.g., by neuromodulators (Burke and Hablitz 1996; Cole and Nicoll 1984; McCormick et al. 1993; Pedarzani and Storm 1993; Spain, 1994)]. Thus the model described here both accounts for the low firing rates observed in networks throughout the mammalian CNS and provides a natural mechanism for switching between steady firing and bursting. A switch between these two patterns is necessary for behavior that is activated episodically, such as locomotion, scratching, and swallowing (Berkinblit et al. 1978; Kudo and Yamada 1987; Zoungrana et al. 1997).

Although endogenously active cells may account for the observed firing patterns in mammalian neuronal networks, they are not necessarily responsible for those firing patterns. This is because the brain is neither isolated nor comprised of a single network. Instead, it receives sensory input and consists of distinct areas that interact through long-range connections. Because external input to a network can drive cells even when they are not receiving input from within the network (making those cells effectively endogenously active), it is possible that the low firing rates observed in cortex are generated by sensory input, or, as has been proposed by Amit and Brunel (1997), by input from other brain areas. The source of the input does not affect our model, however; the model applies whether the input is external or endogenous.

The theory developed here provides a framework for understanding intrinsic firing patterns and their dynamics in large networks of neurons. This does more than simply explain observed firing patterns; it provides a link between basic properties of networks and a way to understand how firing patterns could be modified by patterned external input. This link is critical for developing models of how computations are performed by real neuronal networks.