## Abstract

We explore the mechanism of synchronized bursting activity with frequency of ∼10 Hz that appears in cortical tissues at low extracellular magnesium concentration [Mg^{2+}]_{o}. We hypothesize that this activity is persistent, namely coexists with the quiescent state and depends on slow *N*-methyl-d-aspartate (NMDA) conductances. To explore this hypothesis, we construct and investigate a conductance-based model of excitatory cortical networks. Population bursting activity can persist for physiological values of the NMDA decay time constant (∼100 ms). Neurons are synchronized at the time scale of bursts but not of single spikes. A reduced model of a cell coupled to itself can encompass most of this highly synchronized network behavior and is analyzed using the fast-slow method. Synchronized bursts appear for intermediate values of the NMDA conductance *g*_{NMDA} if NMDA conductances are not too fast. Regular spiking activity appears for larger *g*_{NMDA}. If the single cell is a conditional burster, persistent synchronized bursts become more robust. Weakly synchronized states appear for zero AMPA conductance *g*_{AMPA}. Enhancing *g*_{AMPA} increases both synchrony and the number of spikes within bursts and decreases the bursting frequency. Too strong *g*_{AMPA}, however, prevents the activity because it enhances neuronal intrinsic adaptation. When [Mg^{2+}]_{o} is increased, higher *g*_{NMDA} values are needed to maintain bursting activity. Bursting frequency decreases with [Mg^{2+}]_{o}, and the network is silent with physiological [Mg^{2+}]_{o}. Inhibition weakly decreases the bursting frequency if inhibitory cells receive enough NMDA-mediated excitation. This study explains the importance of conditional bursters in layer V in supporting epileptiform activity at low [Mg^{2+}]_{o}.

## INTRODUCTION

A persistent state of a neuronal network is a state in which periods of increased neuronal discharge dynamically coexist with silent, or an almost silent, epochs. Such neural activity has been demonstrated in the cortex during “working memory” tasks (Funahashi et al. 1989); reviewed in (McCormick et al. 2003). A special type of persistent state, asynchronous firing, was studied using theoretical and computational methods (Amit and Brunel 1997; Hansel and Mato 2001; Wang 1999, 2001). Experiments in monkeys, however, have revealed partially synchronized activity during working memory tasks that have oscillatory (Cardoso de Oliveira et al. 2001) or nonoscillatory (Vaadia et al. 1995) characteristics. The mechanisms that generate persistent synchronized states in cortical networks are still unknown.

Here we study the mechanism for a persistent synchronized state in a specific neocortical system that was studied in vitro: discharges of rhythmic bursting activity at low extracellular Mg^{2+} concentration ([Mg^{2+}]_{o}) solution (Kawaguchi 2001; Sutor and Hablitz 1989). These synchronized oscillations depend on NMDA receptors and appear in tissue segments containing layer V only (Silva et al. 1991). In rats, this layer contains, among other cell types, pyramidal neurons that are quiescent at rest and burst rhythmically when constant, depolarizing current is applied, suggesting that this intrinsic bursting property could facilitate the generation of these waves (Steriade 2004). The oscillation frequency in the somatosensory cortex is ∼8–12 Hz. Blocking GABA_{A}-mediated inhibition considerably increases the amplitude of the local field potentials and the duration of the oscillatory episode (up to ∼3 s) and only slightly increases the oscillation frequency, namely the peak of the power spectrum increases by ∼10–20% (Flint and Connors 1996). With inhibition intact, some inhibitory neurons burst in synchrony with the excitatory cells, whereas other fire almost continuously (Kawaguchi 2001). Recently, properties of these waves were studied using a combination of optical, voltage-sensitive dye imaging and local field potential techniques. Episodes of bursting oscillations can be evoked or can start spontaneously at a few confined initiation foci (Tsau et al. 1998). Correlation between local field potentials measured at two points did not decrease considerably with distance, but there were large phase shifts at a large distance (Fig. 5 in Wu et al. 1999). The phase relationship between the starting times of bursting period was not constant in time (Figs. 4 and 7 in Wu et al. 1999), hinting that the activity after the steady-state activity is a result of collective dynamics and not of periodic stimulation of the slice by a group of pacemakers located at a certain location.

In this paper, we investigate the simplest scenario of evoked propagating waves in a slice model without noise or disorder. We focus on dynamics of order of a few hundred milliseconds or less and do not consider the slow processes of the order of seconds or more that eventually terminate the episode of activity (therefore in reality, the activity is not strictly persistent at long times). First, we study disinhibited networks because inhibitory neurons do not have a major effect on the basic discharge properties (Flint and Connors 1996) of the population spiking activity. Then we explore the effect of inhibition. We hypothesize that the synchronized waves at low [Mg^{2+}]_{o} appear as synchronized persistent bursting activity. This hypothesis raises several questions.

First, the bursting time period of the population (field) is of the order of 100 ms. Because in each burst the active (spiking) phase of the burst eventually terminates, there is a process that terminates it that is slow in comparison with the inter-spike interval within the burst. To generate a subsequent burst, the slow *N*-methyl-d-aspartate (NMDA) conductance should overcome this slow process and initiate the spiking activity. For that, the NMDA-mediated conductance should decay slowly enough. How small can this decay time constant be to still support the next active phase? Is this value within the physiological range?

Second, should the single cells have a propensity to burst to obtain network bursting or can bursting be obtained as a network effect with the single cells only being regular spiking cells?

Third, does the fast, AMPA-mediated excitation support or disrupt persistent bursting? How does it affect the properties of this activity, such as bursting duration, frequency and level of synchronization?

Fourth, how large should [Mg^{2+}]_{o} be to prevent the persistent activity? How do the burst properties depend on [Mg^{2+}]_{o} in the parameter regime where they exist?

And fifth, how does intracortical inhibition affect the discharge frequency? Does the inhibitory effect depend on the type and conductance strength of excitatory receptors (AMPA and NMDA) on inhibitory cells?

To answer these questions, we construct and analyze a conductance-based model of one-dimensional networks with spatially decaying connectivity.

## METHODS

### Single-cell model

##### EXCITATORY CELLS.

We model excitatory cortical cells that can fire either in trains of single action potentials (regular spiking, RS) or in fast bursts of action potentials (intrinsically bursting, IB) (Connors et al. 1982; McCormick et al. 1985). To evaluate the importance of the intrinsic bursting property, we use a model that can be transformed, by modifying a single parameter, from RS to IB with an increasing number of spikes within a burst (Golomb and Amitai 1997; D. Golomb and Y. Yaari, unpublished data). The single-cell dynamics is described by a single-compartment Hodgkin-Huxley type model with the use of a set of coupled differential equations (1) Where *C* is the membrane capacitance and *V*(*x,t*) is the membrane potential of a neuron at a position *x* and time *t*. The right-hand side incorporates an applied current *I*_{app} and the following intrinsic and synaptic currents: the transient sodium current *I*_{Na}, the persistent sodium current *I*_{NaP}, the delayed-rectifier potassium current *I*_{Kdr}, the slow potassium current *I*_{K-slow}, the leak current *I*_{L}, the AMPA current *I*_{AMPA}, the NMDA current *I*_{NMDA}, and the GABA_{A} current *I*_{GABAA}^{IE}. The currents *I*_{Na} and *I*_{Kdr} are needed for spike generation. The slow potassium current *I*_{K-slow} in our model represents potassium currents with kinetics slower than the action potential time scale, i.e., with activation time scale of several tens to hundreds of milliseconds (Storm 1990). These currents are either calcium dependent, such as *I*_{AHP} (Pinsky and Rinzel 1994; Sah 1996), or voltage dependent, such as *I*_{M} (Bibbig et al. 2001; Mainen and Sejnowski 1996; Yue and Yaari 2004). The current *I*_{K-slow} contributes the slow variable that is needed for adaptation or bursting in RS cortical cells (Rinzel and Ermentrout 1998). If the kinetics of *I*_{Na} and *I*_{Kdr} are fast enough, and with strong enough persistent sodium conductance (*g*_{NaP}), the model may produce IB behavior in response to prolonged current injection. For example, for , the neuron bursts if *g*_{NaP} > 0.091 mS/cm^{2}. Neurons from this type are called “conditional intrinsically bursting” (cIB). Hence, the conductance *g*_{NaP} can be viewed as an “intrinsic bursting parameter.” The equations and parameters of the model are given in appendix a.

This single-cell model for cortical pyramidal neurons was chosen for two reasons. First, it is a minimal model that can show both regular spiking and bursting behavior. Second, the role of *g*_{NaP} in the intrinsic bursting property in the model is consistent with experimental observation showing the importance of this conductance for both single-cell bursting in cortical structures (e.g., Su et al. 2001) and for bursting in neocortical networks bathed in low [Mg^{2+}]_{o} (Castro-Alamancos and Rigas 2004).

##### INHIBITORY CELLS.

We model inhibitory cortical cells using the Wang-Buzsáki model of fast spiking (FS) neurons (Wang and Buzsáki 1996). The single-cell dynamics is described by a single-compartment model (2) The equations and parameters of the model are given in appendix a.

### Synaptic models

Neocortical cells receive fast, AMPA-mediated, and slow, NMDA-mediated excitatory postsynaptic potentials (EPSPs) from neighboring excitatory cells (Gil and Amitai 1996), and GABA_{A}-mediated inhibitory postsynaptic potentials (IPSPs) from neighboring inhibitory cells (Hansel and Mato 2003; Wang and Buzsáki 1996). A gating variable *s*_{AMPA} for an AMPA receptor, representing the fraction of open channels, is modeled according to (3) where *V*_{pre} is the presynaptic potential, and (Golomb and Amitai 1997; Wang and Rinzel 1993), . The rise time of the AMPA synapses is assumed to be instantaneous, and decay time of those synapse is (Stern et al. 1992).

The rise time of NMDA receptors cannot be neglected, and therefore two differential equations are needed to model those synapses. We use the phenomenological equations of Golomb et al. (1996) (4) (5) where , , . The decay time of NMDA-mediated excitatory postsynaptic conductances (EPSCs) in cortex varies from 65 to 108 ms (Kumar and Huguenard 2003). For our reference parameter set, we use the value τ_{NMDA} =100 ms, corresponding to the value of intracortically evoked NMDA-mediated EPSCs (Kumar and Huguenard 2003). The effect of this parameter is assessed in the following text. For simplicity, effects of synaptic depression and facilitation (Golomb and Amitai 1997; Tsodyks and Markram 1997) are not included in the model.

A gating variable *s*_{GABAA} for a GABA_{A} receptor, is modeled according to (6) where , (Hansel and Mato 2003; Wang and Buzsaki 1996).

The total synaptic conductance a neuron receives, *gS*^{in} (for AMPA, NMDA, or GABA_{A}) is calculated by summing the synaptic variable of each of its presynaptic neurons as described in the following text. The AMPA current is (7) where *V*_{Glu} is the reversal potential of glutamatergic synapses. If [Mg^{2+}]_{o} is larger than 0, the NMDA current also depends on the postsynaptic voltage and is calculated by multiplying the presynaptic term by a sigmoid function *f*_{NMDA}(*V*) of the postsynaptic voltage (Destexhe et al. 1994; Traub et al. 1991). The NMDA current is therefore (8) The half-maximum voltage of *f*_{NMDA}(*V*) is θ_{NMDA}. The value of θ_{NMDA} is −∞ for and is −31 mV for NMDA-mediated intracortical synapses at physiological levels of [Mg^{2+}]_{o} (Kumar and Huguenard 2003); it increases logarithmically with [Mg^{2+}]_{o} (Jahr and Stevens 1990). Using the data for intracortical connections, we find that θ_{NMDA} = 10.5 mV × ln([Mg^{2+}]_{o}/38.3 mM).

GABA_{A} current is (9) where *V*_{GABAA} is the reversal potential of GABA_{A} synapses.

Axonal propagation delay and synaptic delays are neglected because the time scale of bursting oscillations is in the order of 100 ms, much larger than the time scale of those delays, which are in the order of a few milliseconds.

### Network architecture

Propagation of activity along cortical slices has been modeled before in one-dimensional networks (Ermentrout 1998; Golomb 1998; Golomb and Amitai 1997; Golomb and Ermentrout 1999, 2002; Golomb et al. 2001), where the extensive intralaminar connectivity is neglected, assuming cortical columns are rapidly and fully recruited. In the same manner, our model consists of a one-dimensional system. Excitatory and inhibitory neurons are equally distributed along the interval 0 ≤ *x* ≤ *L*, where *L* is the slice length and *x* is the neuron position. The cell density is ρ, and the number of neurons from each type is . The position of the *i*th neuron, either excitatory or inhibitory, 1 ≤ *i* ≤ *N*, is . The interaction between neurons is assumed to decay with the distance between them. The “synaptic footprint shape,” *w*(*x*), denotes the functional dependence of the synaptic connectivity on the distance between the preand postsynaptic cells. It is assumed here to be exponential with a characteristic delay length (“footprint length”) σ, namely . The model is studied in the parameter regime 1/ρ ≪ σ ≪ *L*. The total synaptic conductance *gS*^{in} affecting a neuron is (10) The quantities *g*, *s*, and *S*^{in} have a subscript AMPA, NMDA, or GABA_{A} for the three types of receptors. They also have superscripts αβ, where α and β and denote the pre- and postsynaptic populations, respectively, that can be E (excitatory) or I (inhibitory). The coupling among neurons in the inhibitory population, either excitatory or inhibitory (Beierlein et al. 2003), is neglected here. To shorten the notation, we omit the superscripts if they are EE. We use open boundary conditions. The integro-differential equations (1–10), together with the differential equations for the auxiliary variables, are discretized on a grid, as described in appendix a.

### Initial conditions

As in our previous work (Golomb 1998; Golomb and Amitai 1997), we initiate our model from a state at which all the neurons, except a group at the left edge, are in their resting state. A wave is initiated by depolarizing a group of neuron, both excitatory and inhibitory, within a length σ or σ/2 at the left edge, such that they generate action potentials. The firing neurons may recruit resting neurons through their synaptic connections to initiate a propagating discharge.

### Definitions of states

We define the state of the system based on its long-time behavior. Therefore a state of transient activity, in which the network is silent after brief transient activity, is considered to be quiescent. A transient bursting state followed by a tonic state is considered to be tonic.

For a single-cell model, or a state of a disinhibited network that can be represented by a single cell coupled to itself, the notion of “bursting cell” means here “a spiking cell that is not firing periodically (tonically).”^{1} To determine whether the voltage time course of an isolated neuron is tonic or bursting, we run a simulation for a long time *T _{m}* after a transient time

*T*

_{transient}. We find the minimal inter-spike interval

*T*

_{isi,min}and the maximal inter-spike interval

*T*

_{isi,max}. The cell is considered to be in a tonic mode if

*T*

_{isi,min}/

*T*

_{isi,max}> 0.9. Otherwise, the cell is in a bursting mode.

The situation is more complicated in a network, where neurons can fire in an irregular and aperiodic manner. Therefore a neuron in a network is defined to be “bursting” if there is a clear distinction between brief and prolonged inter-spike interval, namely *T*_{isi,min}/*T*_{isi,max} < 0.33. The properties of the network are determined by considering the middle group of excitatory neurons with length *L*/2 and ignoring the two edges of length *L*/4. A network is considered to be in a tonic mode if all those middle excitatory neurons fire tonically. Borders between regimes of activity are determined using the bisection method (Golomb and Ermentrout 1999).

For both single neurons and networks, we compute the average number over the excitatory population of spikes within a burst, *N*_{s}, and the bursting frequency, *f*. The burst duration *T*_{d} is defined as the average time between the first spike and the last spike in a burst.

### Correlation calculation.

The average membrane potential of a neuron is (11) The cross-correlation between the membrane potentials of excitatory neurons at positions *x* and *x′* is (12) Note that we subtract the product of the average voltages.

### Synchrony measure

For networks with all-to-all coupling, the synchrony measure χ of the excitatory population is defined as the fluctuation of the population-average membrane potential, normalized by the average fluctuations of the membrane potentials of single neurons (Golomb and Rinzel 1993, 1994; Golomb et al. 2001). We extend this definition for networks with spatially decaying synaptic coupling by calculating the population average over excitatory neurons in the neighborhood of a specific neuron in the middle of the chain, as explained in appendix b. The measure χ varies between 0 (asynchronized state) and 1 (fully synchronized state).

### Fast-slow analysis

We used the fast-slow method (Bertram et al. 1995; Hoppensteadt and Izhikevich 1997; Izhikevich 2000; Rinzel and Ermentrout 1998) to study the bursting mechanism of a single excitatory cell coupled to itself. This method has been applied successfully to analyze periodic bursting of various biophysical models of neurons and networks (e.g., Jung et al. 1996; Mandelblat et al. 2001; Tabak et al. 2000). The method makes use of basic dynamical systems theory (Strogatz 1994), which separates the variables of the system into two subsystems, “fast” and “slow,” according to their overall timescales. In our system (the reduced model of an excitatory cell coupled to itself), the variables *V* (the membrane potential), *h* (inactivation of *I*_{Na}), *n* (activation of *I*_{Kdr}), and s_{AMPA} (the fraction of open AMPA channels), with time scales on the order of 1 ms, are considered to be fast. The variables *z* (activation of *I*_{K-slow}) and s_{NMDA} (the fraction of open NMDA channels), with time scales on the order of 100 ms, are considered to be slow.^{2} In the first stage of the analysis, as will be described later, the variable s_{NMDA} is considered to be constant, and the fast subsystem includes only the variable *z*. The bifurcation diagram of the fast subsystem is computed with the slow variable *z* considered as a parameter. Then the dynamics of the slow subsystem are computed using the time-averaged values of the fast subsystem. This is done by plotting two more curves on the bifurcation diagram. One is the activation curve of *I*_{M}, (*Eq. A15*). The second curve is the equivalent voltage during the oscillatory state, *V*_{equiv} (Bertram et al. 1995; Mandelblat et al. 2001), which is computed to determine the slow evolution of *z* in the full model. The value of *V*_{equiv} for a periodic cycle with a specific *z* is the voltage that, for a rest state (fixed point), yielded a value of d*z*/d*t* that is equal to the value of d*z*/d*t* averaged over the cycle. It was defined implicitly by the equation (13) where *V*(*t*) is the voltage as a function of time along the cycle calculated for a specific value of z and *T _{z}* is the time period of the cycle that depends

*z*. If, for a certain value of z,

*V*

_{equiv}(

*z*) >

*z*

_{∞}

^{−1}(

*z*) (where

*z*

_{∞}

^{−1}(

*z*) means the inverse function of z

_{∞}), than z in the full dynamical system (fast + slow together) increases slowly with a rate of order τ

_{z}. If

*V*

_{equiv}(

*z*) <

*z*

_{∞}

^{−1}(

*z*), z decreases slowly.

The type of bursting we study in this paper belongs to the class of “square-wave bursters” (Rinzel and Ermentrout 1998), where the fast subsystem exhibits bistability of a rest state and a periodic firing state. During the rest state, the slow variable *z* decreases slowly, until the rest state does not exist anymore, and the fast subsystem jumps to the periodic firing state. During this state, the slow variable *z* increases slowly, until the firing state does not exist anymore and the fast subsystem jumps back to its rest state. This mechanism generates periodic bursting.

The bursting dynamics in our system is different from the classical square wave bursting mechanism because the slow variable *s*_{NMDA} is not really a constant. We examine the effects of this second slow variable below.

## RESULTS

We start our investigation of the persistent synchronized bursting states by considering disinhibited networks at . In a second step, we study how the properties of the states vary when [Mg^{2+}]_{o} is elevated. Finally, we explore the role of inhibitory interneurons.

### Synchronized persistent bursting for [Mg^{2+}]_{o} = 0

The disinhibited slice model with the reference parameter set exhibits persistent synchronized bursting activity in response to stimulation at the “left.” The initial burst propagates as a traveling pulse at a constant velocity, similar to propagating single bursts in models of slices in which inhibition was reduced (Golomb and Amitai 1997) (Fig. 1*A*, *left*). Subsequent bursts propagate in a less orderly manner, especially near the edges of the slice model, and do not maintain a constant velocity of the left-to-right propagation. At long times after the initiation of the activity, cells burst with a frequency of , and the bursting activity is synchronized, although not fully synchronized (Fig. 1*A*, *right*). Neighboring cells tend to burst at similar times (Fig. 1*B*, *left*), but spikes within bursts are only weakly synchronized, if at all (Fig. 1*B*, *right*). Summarizing, the model with reference parameter set generates persistent synchronized behavior with properties similar to those seen experimentally (Flint and Connors 1996; Silva et al. 1991; Sutor and Hablitz 1989; Tsau et al. 1998; Wu et al. 1999).

To further quantify the synchronization properties of the bursting activity, we compute the auto-correlation *C*_{x,x}(τ) of the membrane potential *V*(*t*) of one neuron (*Eq. 12*) in the middle of the slice (, Fig. 2*A*) and the cross-correlation *C _{x,x′}*(τ) of this neuron with its neighbor at (Fig. 2

*B*). Both the auto- and cross-correlations oscillate on the bursting time scale, reflecting a periodic firing pattern and synchronization at the bursting level. The auto-correlation shows peaks at the spiking time scale for τ near zero (Fig. 2

*A*,

*top right*) but almost no peaks at the spiking time scale (and only one wide peak at the bursting time scale) for τ value around other integer multiplications of the time period (Fig. 2

*A*,

*bottom right*). The cross-correlation shows almost no peaks at the spiking time scale around (Fig. 2

*B*,

*top right*) and no peaks at the spiking time scale near other integer multiplications of

*T*

_{per}(Fig. 2

*B*,

*bottom right*). For all neurons within a footprint length σ around the neuron at , the amplitude of the cross-correlation is almost constant with the distance between the neurons, and the phase of the cross-correlation is zero. We conclude that the neurons are synchronized on the bursting time scale, at least within distances of order σ, but not on the spiking time scale. Note that noise, heterogeneity, and sparseness are expected to cause the cross-correlations to decay with space and time (Golomb 1998).

The persistent activity in the model does not stop. In the experiments, however, bursting activity in disinhibited slices is terminated after ∼3 s (Flint and Connors 1996). This termination can be explained by the effects of slow processes, such as slow inactivation of Na^{+} channels (Fleidervish and Gutnick 1996; Fleidervish et al. 1996), that are not included in the model.

### Dependence of activity regimes on τ_{NMDA} and g_{NMDA}

The bursting activity of each neuron is terminated by the slow potassium current *I*_{K-slow}, with activation kinetics constant . Once this current has decayed, a strong enough slow NMDA-mediated excitation is responsible for the generation of a new burst. Therefore the decay time constant of the NMDA-mediated synapses, τ_{NMDA}, is an important parameter for generating a persistent synchronized bursting state. We expect that this state can appear only if τ_{NMDA} is large enough in comparison with τ_{z}. Furthermore, one may expect that the τ_{NMDA} value that allows persistent bursting should be large enough in comparison with *T*_{per} and that increasing *g*_{NMDA} may allow persistent bursting even for smaller τ_{NMDA}. To examine the roles of τ_{NMDA} and *g*_{NMDA}, we computed the regimes in the τ_{NMDA}–*g*_{NMDA} plane where stimulation of the slice led to one of three states: persistent bursting activity, persistent tonic activity, or nothing but the quiescent state. The calculation was carried out without (Fig. 3*A*) and with (*B*) the inclusion of AMPA-mediated synaptic conductances (*g*_{AMPA} = 0 and 0.08 mS/cm^{2} respectively). For both cases, bursting activity appear for large enough τ_{NMDA} and moderate values of *g*_{NMDA}. Interestingly, bursting activity is found even for τ_{NMDA} values that are smaller than *τ _{z}*, especially for . For larger values of

*g*

_{NMDA}and for τ

_{NMDA}that is not too small, the activity is tonic, whereas for small

*g*

_{NMDA}or small τ

_{NMDA}, the persistent activity does not prevail. The

*g*

_{NMDA}interval where bursting activity exists increases with τ

_{NMDA}and extends toward low-

*g*

_{NMDA}values. Increasing

*g*

_{AMPA}shifts the boundary of the bursting regime toward larger τ

_{NMDA}values. This means that for moderate values of τ

_{NMDA}, increasing the fast, AMPA-mediated excitation may abolish the persistent bursting behavior. This network behavior stems from the fact that increasing

*g*

_{AMPA}increases the number of spikes within a burst (see following text), and therefore increases the activation of

*I*

_{K-slow}that terminates the burst and prevents the next one. Similarly, the transition from the tonic regime to the bursting regime occurs for larger

*g*

_{NMDA}values as

*g*

_{AMPA}increases, because the stronger

*I*

_{K-slow}due to the increased number of spikes terminate the spiking activity and generates a quiescent period before the next firing period begins.

### Reduced model: cell coupled to itself

The large synchronization on the bursting time scale among neighboring neurons in cases such as those shown in Figs. 1 and 2 allows the use of a simple approximation for the large dynamical system: a neuron coupled to itself with the same *g*_{NMDA} and *g*_{AMPA} values as in the network model (van Vreeswijk and Hansel 2001). This reduced model is amenable to mathematical analysis and is much faster to simulate to produce phase diagrams describing the bursting and tonic regimes as functions of various parameters. This simplified model is only an approximation, however, because the spikes within the bursts are not (or only weakly) synchronized. If the bursts themselves are less synchronized in the full network model, this approximation becomes even worse. Therefore we compare between the behavioral regimes of the reduced model (Fig. 3, *C* and *D*) and the full network model (Fig. 3, *A* and *B*). Qualitatively, the behavioral regimes of the two models are similar. Two main differences are found, however. First, for , the bursting regime extends toward low τ_{NMDA} in the full model in comparison with the reduced one. This effect is caused by the appearance of states with low synchronization in the full model as will be shown in the following text. Second, for , there is a “tristable” regime where the bursting, tonic and quiescent states are stable states. Interestingly, the tristable regime has a tongue-like shape. Each tongue corresponds to a certain number of spikes in the burst, as written in Fig. 3*D*. Because the transition between a certain number of spikes to the subsequent one occurs through chaotic behavior (Terman 1992), it is not surprising that the tongue structure is complex. In our network simulations, the state is selected by the initial conditions, and therefore tristability is not seen. Overall, the reduced model can be used to describe the dynamical regimes of the whole network unless the network exhibits bursts that are not highly synchronized.

### Dependence of bursting patterns on τ_{NMDA} and g_{NMDA}

We next examine how the main properties of the wave: the number of spikes within a burst *N*_{s}, the bursting frequency *f*, the burst duration *T*_{d}, and the synchrony measure χ depend on the properties of NMDA synapses: *g*_{NMDA} (Fig. 4, *A* and *B*) and τ_{NMDA} (Fig. 4, *C* and *D*). The examination is carried out for parameters that lead to bursting, without () and with () AMPA excitation. Computation was performed for both the full network model (red) and the reduced model of one cell coupled to itself (black), to further compare their behavior. In the regimes where the two models yield persistent synchronized bursting, the bursting properties *N*_{s}, *f*, and *T*_{d}, are, in general, very similar, showing that the reduced model mimics the full network model well. The number of spikes *N*_{s} is almost independent of *g*_{NMDA} and τ_{NMDA}. The bursting frequency *f* increases with *g*_{NMDA} but only weakly with τ_{NMDA}. To obtain physiological values of *f*, ∼12 Hz (Flint and Connors 1996), *g*_{NMDA} should not be too large. The burst duration *T*_{d} increases weakly with *g*_{NMDA} and is almost independent of τ_{NMDA}. A main determinant of the burst duration is the time constant of the slow potassium current, *τ _{z}* (see the section on the fast-slow analysis).

### States with low synchrony

For , the full model, but not the reduced model, shows bursting activity for low *g*_{NMDA} values (Fig. 4*A*). These bursting states have low synchrony: χ is <0.2 in comparison to values >0.5 that are obtained for larger *g*_{NMDA} values for which the bursts are highly synchronized. We find that for the value , χ increases sharply. For this same *g*_{NMDA}, *f* of the full model shows a small deflection downward, and the reduced model starts to fire bursts instead of being quiescent. Similarly, states with low bursting synchrony are found for , the reference value of *g*_{NMDA} and low τ_{NMDA}, <60 ms (Fig. 4*C*). For , such low-synchrony states do not appear, and all the bursting states are highly synchronized.

To explore the nature of this state further, we present a rastergram obtained in response to stimulation of the slice model with and in Fig. 5*A*. Initially, the pattern of activity (Fig. 5*A*, *left*) resembles that of the reference parameter set (Fig. 1*A, left*). The firing pattern in the steady state, however, is more complicated (Fig. 5*A, right*). The bursting time of neurons is widely distributed as exhibited both in the rastergram and in the voltage time courses of two bursting neurons (Fig. 5*B*). To determine whether the network is asynchronous, we computed the synchrony measure χ for a local population of 2σρ + 1 neurons, as a function of 1/ρ where ρ varies from 64 to 1,024 (Fig. 5*C*). The value of χ remains roughly constant near the value of 0.17, indicating that the state shows a low level of synchronization and is not asynchronized.

### Dependence of activity regimes and bursting patterns on g_{AMPA} and g_{NaP}

The parameters *g*_{NMDA} and τ_{NMDA} control the slow dynamics of the system, on the order of 100 ms. Bursting behavior, however, depends also on the fast dynamics of the system (Bertram et al. 1995; Izhikevich 2000; Rinzel and Ermentrout 1998). In particular, fast, AMPA-mediated excitation was shown to induce bursting activity in networks of excitatory neurons with spike adaptation (van Vreeswijk and Hansel 2001). Even without AMPA excitation, our model neuron can generate endogenous bursting behavior in response to constant applied current, and the parameter *g*_{NaP} represents the intrinsic conditional bursting property of the single neuron. Without *g*_{NaP}, the neuron fires only tonically in response to a step current pulse. It responds by bursting to this current pulse for large enough *g*_{NaP}, and the number of spikes within each burst increases with increasing *g*_{NaP} (Golomb and Yaari, unpublished).

We study how persistent bursting in the network depends on *g*_{NaP} and *g*_{AMPA} by determining the various behavioral regimes as functions of these parameters using the reduced model. The phase diagram is plotted in Fig. 6 for two values of *g*_{NMDA}: the reference value (*A*) and a value that is 50% larger, (*B*). At the lower *g*_{NMDA} value, tonic firing is obtained for small values of *g*_{NaP} and *g*_{AMPA} (but not for ). Bursts are obtained for moderate *g*_{NaP} values and *g*_{AMPA} values that are not too large. The behavior for the higher *g*_{NMDA} value, 0.105 mS/cm^{2} shares some similarity with the behavior for but differs in the following aspects. First, the tonic regime extends toward larger *g*_{NaP} values. Second, the bursting regime extends toward larger *g*_{AMPA} values. Third, bursts are generated even for , though in a restricted regime. Fourth, a tristable regime emerges for moderate values of *g*_{AMPA} and *g*_{NaP}. For both values of *g*_{NMDA}, increasing *g*_{AMPA} or *g*_{NaP} transfers the activity from the bursting to quiescent states. Reducing *g*_{AMPA} may transfer the activity from bursting to tonic state but not to quiescence (except in the small areas near the tongues). Reducing *g*_{NaP}, however, mostly transfers the activity from bursting to quiescent state. The intrinsic conductance *g*_{NaP} and the synaptic conductance *g*_{AMPA} display therefore some similarities, but also some differences, in their effects on the bursting activity.

The dependence of *N*_{s}, *f*, *T*_{d}, and χ on *g*_{AMPA} is shown in Fig. 7*A* for . Increasing *g*_{AMPA} causes *N*_{s} to increase and *f* to decrease. These two phenomena are related because an increase in *N*_{s} leads to an increase in the activation variables of *I*_{K-slow}, resulting in an elongation of the inter-burst interval. The burst duration *T*_{d} increases only weakly with *g*_{AMPA}. Larger *N*_{s} values that result from increasing *g*_{AMPA} indeed prolong the burst duration, but when *N*_{s} is kept constant, the burst duration even decreases with *g*_{AMPA}. The bursts are synchronized (χ > 0.5), and χ gradually increases with *g*_{AMPA} with many local peaks corresponding to minima in *T*_{d}.

The effect of *N*_{s} and *f* on *g*_{NaP} is similar to the effect of *g*_{AMPA}, whether *g*_{AMPA} is zero (Fig. 7*B*) or *g*_{AMPA} = 0.08 mS/cm^{2} (Fig. 7*C*). The frequency *f*, however, depends only weakly on *g*_{NaP} for , and *T*_{d} is almost constant with some fluctuations. For , χ fluctuates and may reach low values representing states of low synchrony (Fig. 7*B*, *bottom*). Such low synchrony states do not exist for larger *g*_{AMPA} (Fig. 7*C*, *bottom*). We conclude that the network exhibits persistent burst synchronization with substantial AMPA excitation. Without it, bursting states of low synchrony appear in many parameter regimes.

### Fast-slow analysis of bursting activity

##### DESCRIPTION OF THE ANALYSIS.

The reduced model of one cell coupled to itself describes the properties of the network discharges very well (Fig. 3, 4, and 7). Therefore we use this simple model to analyze the properties of the persistent synchronized activity, using the fast-slow method (Bertram et al. 1995; Izhikevich 2000; Rinzel and Ermentrout 1998). We separate the variables of the system (the reduced model) into two subsystems, fast and slow, according to their overall time scales. The fast subsystem includes the variables *V*, *h*, *n*, and s_{AMPA}, and the slow subsystem includes the variables *z* and s_{NMDA}. We analyze the behavior of the single neuron coupled to itself for the reference parameter set, *g*_{NMDA} = 0.07 mS/cm^{2}, *g*_{AMPA} = 0.08 mS/cm^{2} (Fig. 8, *A*–*C*). To describe the effect of the two types of synapses on the burst properties, we carry out the analysis also for two more cases: AMPA excitation is blocked (*g*_{NMDA} = 0.07 mS/cm^{2}, *g*_{AMPA} = 0; Fig. 8, *D*–*F*) and NMDA excitation is enhanced (*g*_{NMDA} = 0.105 mS/cm^{2}, *g*_{AMPA} = 0.08 mS/cm^{2}; Fig. 8, *G*–*I*). As a first step, we plot the voltage time course, *V*(*t*), and the time courses of the slow variables *z*(*t*) and *s*_{NMDA}(*t*), computed by simulating the reduced model during one time period *T*_{per}, in Fig. 8, *A*, *D*, and *G*, for those three parameter sets.

When a cell bursts, its state alternates between an active firing state and a silent rest state. The slow hyperpolarizing variable *z* and the slow depolarizing variable *s*_{NMDA} increase during the firing phase of the burst and decrease during the silent phase. Their role in the burst generation, however, is opposite. As *z* increases, it terminates the firing phase of burst, and as it decreases, it may allow firing to resume again. In contrast, *s*_{NMDA} delays the termination of the firing phase. If *s*_{NMDA} decreases too much during the silent phase, the neuron does not resume firing. Furthermore, increasing *τ*_{NMDA} >100 ms affects the bursting properties only weakly, namely the behavior for physiological value of τ_{NMDA} is quite similar to the behavior for very large *τ*_{NMDA}. Therefore as a first step, we study the system in the limit τ_{NMDA}≫ τ_{z} and fix the slow variable *s*_{NMDA} at a constant, average value it obtains during the cycle: 0.85, 0.92, and 0.93 for the three parameter sets (Fig. 8, *A*, *D*, and *G*, respectively, *bottom*) and consider the slow variable *z* to be a parameter. We compute the bifurcation diagrams of the fast subsystem as a function of *z* (Fig. 8, *B*, *E*, and *H*). The steady state (fixed point; thin black line) is stable for large *z*. This stable rest state coalesces with an unstable state and ceases to exist in a saddle-node bifurcation. The rest state, this time a high plateau, is stable again for small values of *z*. At the *z* value where the high rest state gains its stability (a Hopf bifurcation), an oscillatory state (limit cycle) emerges, corresponding to periodic firing. This oscillatory firing state extends toward the right. We plot the minimal and maximal values of the membrane potential during the firing state (thick black lines). There is a narrow regime in which the continuity of the stable oscillatory state is disturbed by an unstable oscillatory state, and for large values of *z* the amplitude of the oscillatory state that gains its stability is larger than the oscillation amplitude for smaller *z*. The oscillatory state coalesces with another unstable oscillator state (saddle-node of periodics bifurcation) that immediately terminates by intersecting the unstable fixed point (homoclinic, or saddle-loop, bifurcation).

Bistability of firing (oscillatory) and quiescent states in the dynamics of the fast subsystem is necessary, but not sufficient, for allowing bursting activity with one slow variable. To further determine the conditions for bursting in the full (fast and slow) model of a cell coupled to itself, we plot two more curves: the activation curve of the current *I*_{K-slow,} *z*_{∞}(*V*) (*Eq. A15*), plotted in green, and the equivalent voltage during the oscillatory state of the fast subsystem, *V*_{equiv} (*Eq. 13*), plotted in blue. For values of *V*_{equiv} above this green line, *z* slowly increases, and for values of *V*_{equiv} below this green line, *z* slowly decreases. For the three parameter sets we investigate (Fig. 8, *B*, *E*, and *H*), the solid (stable) blue line is above the green line, and therefore *z* increases during the firing state until it terminates the oscillations. The thin black line is below the green line, and therefore *z* decreases during the silent state until the rest state disappears (at the saddle-node bifurcation). The full dynamical system exhibits, therefore periodic bursting (“square-wave bursting”) (Rinzel and Ermentrout 1998) for the three parameter sets shown in Fig. 8.

To complete the analysis, the two variables should be considered as slow variables with similar time scale. The activity regimes: rest, bistability, and firing are shown in the *s*_{NMDA}-z plane (Fig. 8, *C*, *F*, and *I*). The boundaries of those regimes are diagonal lines, and their *z* coordinate increases with their *s*_{NMDA} coordinate. The projection of the trajectory of the bursting state is shown in green. It oscillates through the bistable regime between the silent, rest state and the firing state. Both intrinsic neuronal properties (Fig. 8*E*) and AMPA-mediated excitation (Fig. 8*B*) contribute to the generation of this bistable regime. Using this analysis, we now explain the dynamical behavior of the model.

##### EXPLAINING THE EFFECTS OF VARYING *g*_{AMPA}.

We first explain why *N*_{s} increases with *g*_{AMPA} and *f* decreases with it, as shown in Fig. 7*A*. Comparison between the two values of *g*_{AMPA} (Fig. 8, *B* and *E*) shows that the bistable regime is wider for larger *g*_{AMPA}. The low rest states (both stable and unstable) are not affected by *g*_{AMPA} because the AMPA activation curve is substantially larger than 0 only for very depolarized potentials (the dependencies of the rest voltages on *z* are slightly different for the two values of *g*_{AMPA} because of the different values of the constant *s*_{NMDA}). The conductance *g*_{AMPA}, however, affects the periodic action potentials by increasing their frequency and reducing their amplitude. This happens because after the self-coupled neuron fires an action potential, the depolarizing conductance advances the firing of the subsequent one, the sodium current *I*_{Na} is still partially inactivated, and therefore the action potential overshoot is lower, the activation of *I*_{Kdr} is smaller, and the afterhyperpolarization that follows the action potential is less pronounced. Hence, for the same *z*, the minimal voltage during the action potential is more depolarized for larger *g*_{AMPA}. The value of *z*, for which the minimal value of the action potential nearly coincides with the voltage of the unstable rest state, is larger and the bistable regime is wider when *g*_{AMPA} increases. As a result, *N*_{s} is larger and *f* is smaller for larger *g*_{AMPA}.

The intrinsic conductance *g*_{NaP} also increases the bistable regime (not shown), and therefore *N*_{s} increases with *g*_{NaP} (Fig. 7, *B* and *C*).

##### EXPLAINING THE EFFECTS OF VARYING τ_{NMDA}.

We now turn to explain the dependence of *N*_{s} on τ_{NMDA} and why activity stops at small τ_{NMDA} values. The analysis becomes somewhat more complicated by the fact that the value of τ_{NMDA} is not much larger than τ_{z}. If the kinetics of *s*_{NMDA} were very slow, the direction of the projection of the neuronal dynamics onto this plane would be parallel to the ordinate (*z*-axis; Fig. 8*F*). Because these kinetics are not much slower than that of *z*, the value of both *z* and *s*_{NMDA} increase during the firing phase and decrease during the silent phase. Because the boundaries of the bistable regime are diagonal, the *z* interval that the neuronal trajectory covers during the bursting cycle is larger for faster NMDA kinetics. As a result, the number of spikes decreases somewhat with *τ*_{NMDA}, for small *τ*_{NMDA} values, as shown in Fig. 4, *C* and *D*.

If τ_{NMDA} becomes smaller in comparison to τ_{z}, the trajectory of the dynamical system in the *s*_{NMDA}-*z* plane becomes more and more parallel to the abscissa. After firing several spikes, the system is stuck at the rest state (on the left), the bursting process is terminated, *s*_{NMDA} decreases to 0 and *z* decreases to a very small value. This explains why persistent bursts cannot occur for too small τ_{NMDA} (Fig. 3).

##### EXPLAINING THE EFFECTS OF VARYING *g*_{NMDA}.

The numerical simulations (Figs. 4, *A* and *B*, and 3) show that, first, in the persistent bursts regime, *f* increases with *g*_{NMDA}, whereas *N _{s}* and

*T*

_{d}are almost independent of

*g*

_{NMDA}; second, for large enough

*g*

_{NMDA}, the activity becomes tonic instead of bursting. To explain these two numerical observations, we note that the slow currents

*I*

_{K-slow}and

*I*

_{NMDA}for the self-coupled neuron at are (

*Eqs. A13*and

*A20*) (14) (15) If

*z*and

*s*

_{NMDA}are considered constants, varying

*g*

_{K-slow}×

*z*or

*g*

_{NMDA}×

*s*

_{NMDA}is roughly like varying the external current

*I*

_{app}. There is, however, a difference that stems from the driving force term

*V*−

*V*

_{rev}, where

*V*

_{rev}, the reversal potential of K

^{+}or glutamate. Because of this driving force, increasing

*z*is more powerful for depolarizing potentials, whereas increasing

*s*

_{NMDA}is more powerful for hyperpolarizing potentials. This means that when

*g*

_{NMDA}×

*s*

_{NMDA}increases, the bifurcation diagram should move to the right (toward larger

*z*values), and it does it more for hyperpolarizing currents to compensate for the smaller driving force. This is why both the steady-state curve and the firing curves are shifted to the right in Fig. 8

*H*in comparison with Fig. 8

*B*, but the steady-state curve is shifted more to the right and the bistable regime is smaller.

If *g*_{NMDA} is large enough, the solid blue curve, representing the effective potential *V*_{equiv} for the firing state, is shifted enough to the right such that it intersects with the green curve. In this case, the firing state is stable (in the full system) and the system does not burst.

If *g*_{NMDA} is large, but still below the value for the tonic behavior of the full system, the neuron still bursts but the bistable regime is narrower. The value of d*z/*d*t* in the full system depends on [*z*_{∞}(*V*) − *z*)/τ_{z} (*Eq. A14*), namely on the distance between the thin solid black line and the green line for rest states, or the solid blue line and the green line for the tonic states. As shown in Fig. 8*H*, the impact of the shift to the right of the steady state curve is weak, because the value *z*_{∞}(*V*) − *z* remains about 0.1. The bistable regime is narrower, and therefore the silent phase of the burst decreases as *g*_{NMDA} increases. The shift of the solid blue curve to the right is much more important because the solid blue curve and the green curve almost intersect. Therefore the value of d*z*/d*t* during the firing phase of the burst is smaller for larger *g*_{NMDA}. This compensates for the decrease in the *z* regime where the system is bistable, and therefore the burst duration *T*_{d}, as well as *N _{s}*, are almost independent of

*g*

_{NMDA}.

### Effects of nonzero [Mg^{2+}]_{o}

Raising [Mg^{2+}]_{o} leads to several physiological effect. The main effect is on the response of NMDA receptors. Under physiological conditions, when the membrane potential is close to its resting values, NMDA receptors are blocked by Mg^{2+} ions. Increasing [Mg^{2+}]_{o} corresponds to shifting the sigmoid function *f*_{NMDA}(*V*), the voltage dependency of the NMDA synaptic conductance on the postsynaptic voltage, toward the right. This means that θ_{NMDA}, the half-maximum value of *f*_{NMDA}(*V*), is depolarized as [Mg^{2+}]_{o} increases from −∞ for to around −30 mV for physiological value (Jahr and Stevens 1990; Kumar and Huguenard 2003). In addition, elevating [Mg^{2+}]_{o}, like elevating the concentration of other divalent cations, decreases synaptic transmission (Johnston and Wu 1995) and increase the action potential threshold (Frankenhaeuser and Hodgkin 1957). We analyze the consequences of varying θ_{NMDA} and then address the two other effects.

##### EFFECTS OF VARYING θ_{NMDA}.

We computed the dynamical regimes in the θ_{NMDA}-*g*_{NMDA} plane for (Fig. 9*A*) and for (Fig. 9*B*). The bursting regime has a “backward L” shape. Whereas at hyperpolarized θ_{NMDA} values, it gradually shifts to larger *g*_{NMDA} values, as θ_{NMDA} increases, it turns toward very large *g*_{NMDA} values for θ_{NMDA} above −60 mV. The firing patterns at more depolarized θ_{NMDA} and large *g*_{NMDA} are characterized by bursting with a large number of fast spikes in a burst (see traces on the right of Fig. 9, *A* and *B*). No bursting is obtained for *g*_{NMDA} values that are not extremely large and θ_{NMDA} above −40 mV for or −55 mV for . Thus the model predicts that there is no persistent synchronized activity under physiological conditions for (Kumar and Huguenard 2003). Moreover, raising *g*_{AMPA} is expected to shift the bursting regime toward more negative θ_{NMDA} values.

The dependence of *N*_{s}, *f*, *T*_{d}, and χ on θ_{NMDA} is shown in Fig. 9*C* for . The number of spikes *N*_{s} increases and the frequency *f* decreases with θ_{NMDA}. This means that, for nonzero [Mg^{2+}]_{o} values, neurons fire low-frequency bursts with large number of spikes. The burst duration *T*_{d} depends only weakly on θ_{NMDA}, and the system is in a highly synchronized state (χ above 0.5).

##### EXPLAINING THE EFFECTS OF VARYING θ_{NMDA}.

We now explain, using the fast-slow analysis, the dependence of *N*_{s}, *f*, *T*_{d}, and χ. We perform this analysis of the reduced model of one cell coupled to itself for , and (Fig. 10) and compare this parameter set with a similar one with (Fig. 7, *G–I*). Because the frequency *f* is lower for more depolarized θ_{NMDA}, *s*_{NMDA} decreases to lower values (Fig. 10*A*, *bottom*), and hence the bifurcation diagram (Fig. 10*B*) is computed for an average value of . The structure of the bifurcationdiagram for resembles that for θ_{NMDA} = −∞, but there are important quantitative differences. Because raising θ_{NMDA} has a stronger impact on the NMDA current at hyperpolarized voltages (the effect of Mg^{2+} blockade on NMDA-mediated excitation is more effective at hyperpolarized potentials), the bifurcation diagram is shifted more to the left at more negative values of *V*. This means that both the rest state and the firing state extend to lower values of *z*, but the rest state extends more, and therefore the bistable regime is wider for more depolarized θ_{NMDA}. In addition, the rest state potential is more hyperpolarized because the depolarizing effect of the NMDA excitation is weaker because of the Mg^{2+} blockade. The widening of the bistable state and its shift toward lower values of *z* are not just a consequence of the smaller value of *s*_{NMDA} for more depolarized θ_{NMDA}, as shown in Fig. 10*C* (compare with Fig. 7*H*). Because the bistable regime is wider and the rest state is more hyperpolarized, the silent phase of the burst is more prolonged and reaches lower values of *V*. Hence, *f* decreases with θ_{NMDA}. The neuron fires more spikes with θ_{NMDA} because the bistable regime extends toward the *z* values for which the firing state has lower amplitude and higher firing frequencies. In contrast to the duration of the silent phase of the burst, the duration of the firing state *T*_{d} remains more or less the same despite the widening of the bistable regime. As in the case of increasing *g*_{NMDA}, this is a result of the fact that shifting the firing state to the right increases the distance between the blue and green lines, and therefore the value of d*z/*d*t* during the firing phase (*Eq. 13*) is larger for θ_{NMDA} = −70 (Fig. 10*B*) than for (Fig. 7*H*). The dynamics of *z* during the firing phase of the burst is therefore faster and compensates for the widening of the bistable regime, such that *T*_{d} remains almost constant.

##### DIVALENT CATION EFFECTS.

Elevating [Mg^{2+}]_{o} increases firing threshold. We mimic this effect phenomenologically in the model by decreasing the rest potential of the neuron, namely by hyperpolarizing *V*_{L}. This causes a small increase in *N*_{s} and a decrease in *f*; *T*_{d}, and χ are almost unchanged (Fig. 11). If *V*_{L} is too depolarized, the persistent activity ceases to exist. Hence, hyperpolarizing *V*_{L} affects the activity in a similar way to depolarizing θ_{NMDA}.

Decreasing synaptic transmission as a result of elevating [Mg^{2+}]_{o} leads to a decrease in *g*_{AMPA} and *g*_{NMDA}. Decreasing *g*_{NMDA} causes *f* to decrease, almost does not affect *N*_{s}, and eventually terminates the activity (Fig. 4, *A* and *B*). Its effect is therefore similar to the effect of depolarizing θ_{NMDA}. Decreasing *g*_{AMPA} has an opposite effect: it causes *f* to increase and *N*_{s} to decrease (Fig. 7, *A* and *B*) and extends the regime in which persistent bursting exist toward larger θ_{NMDA} values (Fig. 9*A*). Note that persist bursting activity cannot occur for θ_{NMDA} = −31 mV even for *g*_{AMPA} = 0. Near the parameter regime in which bursting disappear, the bursting frequency is large even for .

### Networks with inhibition

Most of the experimental recordings in cortical slices with low [Mg^{2+}]_{o} were carried out with inhibition intact. We therefore examine the effect of including an inhibitory population on the dynamics of the persistent synchronized bursting. An example is shown in Fig. 12, *A* and *B*, where the excitatory cells with the reference parameter set (Fig. 1) are inhibited by GABA_{A}-mediated receptors with conductance strength . The inhibitory cells are excited by AMPA-mediated receptors with . Rastergrams of the excitatory and inhibitory neurons are shown in Fig. 12*A*, and voltage time traces of one excitatory neuron and one inhibitory neuron in the middle of the slice model are shown in Fig. 12*B*. The firing pattern of the excitatory cells is similar to the firing pattern without inhibition (compare Fig. 12, *A* and *B*, with Fig. 1). The inhibitory cells burst as well, and their bursting times follow those of the adjacent excitatory cells. The firing patterns of the inhibitory cells, but not those of excitatory cells, are different when the inhibitory cells are also excited by NMDA-mediated receptors with = 0.05 mS/cm^{2} (Fig. 12, *C* and *D*). The sustained NMDA-mediated excitation causes the inhibitory cell to fire even during the silent periods of the excitatory cells, although the firing rate of the inhibitory cells is enhanced during the firing period of the excitatory cells in response to the fast AMPA-mediated excitation. The number of spikes *N _{s}* fired by excitatory cells during a burst decreases with increasing the excitatory-to-inhibitory synaptic conductances (Fig. 13,

*A*and

*B*) or the inhibitory-to-excitatory synaptic conductances (

*C*and

*D*) because the recurrent inhibition contributes to the burst termination. This decrease in

*N*causes an increase in

_{s}*f*when because the slow K

^{+}current is less activated (Fig. 13,

*A*and

*C*). When

*g*

_{NMDA}

^{EI}is large enough and the inhibitory neurons fire during the whole cycle, they hyperpolarize the excitatory cells and reduce

*f*(Fig. 13,

*B*and

*D*). To conclude, our modeling results show that GABA

_{A}-mediated inhibition causes the

*f*to increases gradually if the inhibitory neurons are excited by AMPA-mediated synaptic conductance only and to decrease gradually if the inhibitory neurons are also excited by strong enough NMDA-mediated synaptic conductance.

## DISCUSSION

We have shown that a network of excitatory neurons coupled by NMDA-mediated synapses can show persistent synchronized bursting activity. Our model proposes two main physiologically relevant conditions for the emergence of such activity: *1*) the presence of enough neurons with intrinsic bursting properties in the population, or enough AMPA excitation; and *2*) relatively open NMDA receptors, such that may exist if θ_{NMDA} is hyperpolarized or the neuronal population is depolarized enough. This persistent activity is robust for physiological values of the τ_{NMDA} (decay time constant of NMDA-mediated synapses) of ∼100 ms. In addition, the conductance *g*_{NMDA} should be moderate; large *g*_{NMDA} supports tonic firing. AMPA conductance often shifts tonic to bursting activity. Surprisingly, *g*_{AMPA} that is too large abolishes all persistent activity in the model. Bursts are often highly synchronized, whereas the spikes within the burst are not (Pinsky and Rinzel 1994; van Vreeswijk and Hansel 2001). Without AMPA conductance, the system may settle into a weakly synchronized state. Elevating [Mg^{2+}]_{o} decreases the bursting frequency and eventually terminates the activity. Even strong increases in *g*_{NMDA} cannot generate synchronized persistent activity in the cortical slice for physiological [Mg^{2+}]_{o}. If inhibitory cells do not receive NMDA-mediated excitation, inhibition increases *f*; if they receive strong NMDA-mediated excitation, inhibition decreases *f*.

### Consequences of the specific model

We have analyzed a bursting neuron model from the “square-wave” type (Bertram et al. 1995; Izhikevich 2000; Rinzel and Ermentrout 1998), which is based on bistability of the fast subsystem of variables. By applying the fast-slow analysis, we have shown that the results are a consequence of the bursting mechanism and not just of the specific choice of intrinsic and synaptic parameters. We expect most of the main qualitative results of this article pertaining to the persistent synchronized bursts to be valid also for other bursting mechanisms. For every burster, τ_{NMDA} should be large enough to overcome the hyperpolarizing current(s) that terminate a burst and start a new one. The exact minimal value of τ_{NMDA} is parameter-dependent, but roughly it should be at least around the time period, 100 ms. This value is similar to the decay rate of NMDA-mediated EPSCs in intracortical synapses (Kumar and Huguenard 2003). If *g*_{NMDA} is too large, the activity is not terminated and neurons fire tonically. Increasing AMPA is expected to increase the number of spikes within a burst and therefore terminate the activity completely in models where a burst is terminated by activation of hyperpolarizing currents.

One necessary condition for this model to generate bursts, however, may not be needed in other bursting mechanisms. In our model, the kinetics of the spike-generating currents *I*_{Na} and *I*_{Kdr} are fast enough to allow bistability in the model of a cell coupled to itself (Bertram et al. 1995; Rinzel and Ermentrout 1998). In other bursting models, for example, those currents based on ping-ponging between soma and dendrite (Mainen and Sejnowski 1996; Pinsky and Rinzel 1994) may have slower kinetics.

### Roles of g_{AMPA} and g_{NaP}

The parameter *g*_{NaP} controls the tendency of the isolated neuron to burst in response to external depolarizing current. According to our model, for small or moderate values of *g*_{AMPA} and moderate values of *g*_{NMDA}, moderate values of *g*_{NaP} are needed to obtain network bursting. This result is consistent with the experimental observation that drugs that suppress *I*_{NaP} abolish the 10-Hz oscillations in the agranular neocortex (Castro-Alamancos and Rigas 2004).

The AMPA conductance *g*_{AMPA} may support persistent bursting behavior even when , but very large *g*_{AMPA} values terminate the activity (Fig. 6). This conductance shifts the bursting regime toward larger *g*_{NMDA} values (Fig. 3). In addition, enhancing *g*_{AMPA} increases the number of spikes within a burst, decreases the bursting frequency, enhances bursting synchronization, and prevents the appearance of weakly synchronized bursting states (Fig. 7).

In the reduced model of a cell coupled to itself, the synaptic AMPA conductance and the persistent Na^{+} (NaP) conductance have a similar, although somewhat different, properties. Their kinetics are fast, but that of the NaP conductance is much faster and is considered here to be instantaneous. They are activated at depolarized membrane potential, but the NaP conductance is activated above , whereas the AMPA conductance is activated above (namely, when the presynaptic neuron fires an action potential). As a result, the dynamical roles of the two conductances, while similar, have some differences. Increasing *g*_{NaP}, like increasing *g*_{AMPA}, extends the bistable regime in the bifurcation diagram of the fast subsystem and therefore causes an increases in the number of spikes within a burst and, in general, a reduction of bursting frequency. Because the NaP conductances are activated above −47 mV, increasing *g*_{NaP}, but not *g*_{AMPA}, can support a high plateau or fast spiking riding on a high plateau. Decreasing *g*_{NaP}, but not *g*_{AMPA}, can transfer a bursting state into a quiescent state (Fig. 6). This is a result of the fact that the NaP conductance is partially activated in subthreshold voltages and can support the firing of a new spike, whereas the AMPA conductance is activated much above firing threshold.

The reduced model cannot, of course, address the issue of synchronization in large networks. Indeed, the AMPA conductance supports synchronization and prevents weakly synchronized states, whereas such states may exist with substantial *g*_{NaP} values.

### Contribution to the theory of synchronization and comparison with other models

Most theoretical works on synchronization deal with synchrony of spikes (reviewed in Golomb et al. 2001). Only a few models deal with synchronization of bursts of spikes, and in most of them, the activity is not persistent, namely there is no stable rest state that coexists with the firing state. Our model differs from those models by creating an active state that coexists with a quiescent state. Still, it is possible to compare the bursting mechanism of our model with the bursting mechanism in those models by noticing that the applied current there plays a similar role to the role of NMDA synaptic current here, which is roughly constant over the duration of a burst cycle. Previous work (van Vreeswijk and Hansel 2001) describes bursting activity that emerges from an excitatory network of neurons as a result of adaptation in neurons that are firing tonically when isolated. Similarly to our model, bistability of the fast subsystem (without the activation variable of the adaptation current, which is similar to *z*) is needed for generating bursting activity. Such bistability is generated in the model of van Vreeswijk and Hansel (2001) by the effect *g*_{AMPA} only, whereas in our model, it can be generated by *g*_{AMPA} or *g*_{NaP}. The active phase of the burst in the conductance-based model of van Vreeswijk and Hansel (2001) is characterized by the strong undershoot of the membrane potential, which is more hyperpolarized than the membrane potential during the silent phase of the burst (Fig. 10*B* there). In our model, the minimal membrane potential during the active phase of the burst is more depolarized than the membrane potential during the silent phase of the burst (Figs. 1, 5, 8, and 9) as found in intracellular recording from cortical neurons during the NMDA-induced activity (Kawaguchi 2001; Sutor and Hablitz 1989).

Like our model and the model of van Vreeswijk and Hansel (2001), highly, but not fully, synchronized bursts can be obtained in networks of spontaneous bursters (Pinsky and Rinzel 1994). In addition, we find states with low, but not zero, synchrony that occur with slow excitation only.

A model of epileptiform activity induced by low Mg^{2+} in the rat hippocampal slices (Traub et al. 1994) showed bursting activity for hundreds of milliseconds. NMDA receptors in that model decayed very slowly, with the time scale of desensitization (Fig. 5 in Traub et al. 1994). In our work, we study the roles of NMDA kinetics, synaptic and intrinsic strengths in determining bursting mechanism, and level of synchrony, issues that were not addressed there. Persistent tonic activity was found in a model of NMDA-coupled neurons without Mg^{2+} blockade (Ermentrout 2003). That model does not show bursting because the neurons lacked AMPA excitation and the necessary intrinsic properties, for example, slow adaptation.

Bursting activity can be generated in networks of excitatory and inhibitory neurons without adaptation (Hansel and Mato 2003) as a result of destabilization of an asynchronous firing state. Those synchronized bursts are not persistent. It can be shown, however, that the persistent synchronized bursting states can occur in the model of Hansel and Mato (2003) if excitatory synapse are slow enough (D. Hansel, private communication). The inhibitory population in that network model takes the role of adaptation in our model in terminating the firing phase of the burst. Unlike our model, such a model will not be consistent with the experimental observation that persistent synchronized states occur in disinhibited networks.

Slow NMDA-mediated synapses without Mg^{2+} blockade were shown to be important for generating persistent asynchronized activity in excitatory networks (Compte et al. 2000; Wang 1999) but not necessarily when an inhibitory population is added to the model (Hansel and Mato 2001). Here we have shown that when neurons possess an adaptation current and the have the proper intrinsic and/or synaptic properties (namely, bistability of the fast subsystem), NMDA-mediated synapses can lead to synchronized bursting activity. Mg^{2+} blockade of NMDA synapses prevents the appearance of persistent bursting for physiological [Mg^{2+}]_{o}. It will be interesting to systematically examine the effects of [Mg^{2+}]_{o} on persistent asynchronized states.

Endogenously active cells may spontaneously initiate a bursting episode (Latham et al. 2000). When a neuronal system shows bistability of a bursting state and a quiescent state, as described here, an episode of bursting activity can also be evoked by an external stimulus. Bursting activity, without bistability with a quiescent state, is observed in cultures of cortical neurons receiving stochastic input (Giugliano et al. 2004).

### Comparison with experimental results.

The bursting frequency *f* in disinhibited slices is ∼10–14 Hz (Flint and Connors 1996). Our model predicts that *f* cannot be much less than that because in that case the NMDA conductance would decay and nothing would start a new burst after the previous one has terminated. The results of our model (Figs. 3 and 4) show that decreasing *f* (or, equivalently, decreasing τ_{NMDA}) cannot be compensated by increasing *g*_{NMDA} because tonic, and not bursting, activity is obtained for large *g*_{NMDA}. With AMPA-mediated synapses blocked and GABA_{A} inhibition intact, Flint and Connors (1996) recorded only two transient population bursts within a time difference of ∼300 ms. This activity is classified as transient and not as persistent according to the classification of this paper because it does not last for many cycles.

Local field potential recording and optical imaging reveal that the bursting activity is synchronized locally (Silva et al. 1991) and that phase shifts can develop as the distance between neurons grows (Wu et al. 1999). Similar behavior is obtained in our model (Figs. 1 and 2). The phase shifts are expected to grow if the network is heterogeneous. The fact that the activity in our model is persistent and does not result from entrainment by a group of pacemakers is consistent with the optical imaging results that the phase relationship between starting times of bursts vary with time (Wu et al. 1999).

Our modeling results indicate that inhibitory neurons fire bursts in synchrony with their neighboring excitatory cells if they receive mainly AMPA-mediated excitation; they fire almost continuously if they receive strong enough NMDA-mediated excitation. Intracellular recordings from fast-spiking inhibitory interneurons show examples for the two firing patterns (Figs. 4*D* and 8*A* in Kawaguchi 2001). Even inhibitory cells that fire continuously tend to fire more spikes during the bursting period of the excitatory cells (Fig. 8 in Kawaguchi 2001), in agreement with our simulations (Fig. 12). Our theory suggests that blocking GABA_{A}-mediated inhibition decreases the frequency *f* of the network bursting if inhibitory neurons are excited mainly by AMPA-mediated receptors, but this blockade increases *f* if inhibitory neurons are substantially excited by NMDA-mediated receptors (Fig. 13). Because experiments carried out in the somatosensory cortex revealed that the bursting frequency increases weakly with blocking inhibition (Flint and Connors 1996), we predict that inhibitory cells in this cortical area should receive NMDA-mediated excitation from neighboring excitatory cells.

Our model takes into account kinetic processes that are faster than ∼100 ms. Slow processes, such as slow inactivation of Na^{+} currents (Fleidervish and Gutnick 1996; Fleidervish et al. 1996) or slow synaptic depression, are neglected. As a result, the synchronized bursting activity is persistent and does not stop. In the experimental systems, those slow processes and others finally terminate the activity after <1 s if inhibition is intact and after ∼3 s if inhibition is blocked (Flint and Connors 1996). The epileptic-like activity generated in response to lowering [Mg^{2+}]_{o} leads also to accumulation of extracellular Ca^{2+} and K^{+} ion and the alternation of the K^{+}-Cl^{−} electrogenic pump. Such slow processes are not considered here.

Short-term synaptic processes, depression and facilitation, are neglected here for simplicity (Golomb and Amitai 1997; Markram and Tsodyks 1996). Because NMDA-mediated receptors in this model are close to saturation after a brief burst of a few spikes, the effects of depression and facilitation on the NMDA-mediated excitation are expected to be weak. Depression effects on AMPA-mediated excitation can be mostly compensated for by increasing *g*_{AMPA}. The residual effect may lead to a moderate decrease in *N*_{s}.

The present paper describes a deterministic and ordered network. Noise, heterogeneity, and sparseness are expected to cause the cross-correlations to decay with space and time and to eliminate any synchrony at the spiking time scale (Fig. 2). Only large levels of sparseness are expected to desynchronize bursts of neighboring neurons (Golomb 1998).

The most pronounced effect of elevating [Mg^{2+}]_{o} is a shift of the activation curve *f*_{NMDA}(V) (*Eq. 8*) toward more depolarized values, namely making θ_{NMDA} more depolarized. As a result, *f* decrease and *N*_{s} increases, and for large enough θ_{NMDA} the persistent activity cannot exist (Fig. 9). Elevating [Mg^{2+}]_{o} has several other effects, although the strength of these effects in cortical neurons is not known quantitatively. It causes an increase in firing threshold and decrease of synaptic transmission. The effects of raising firing threshold and reducing NMDA-mediated synaptic conductance enhance the effects of depolarizing θ_{NMDA}. Reducing AMPA-mediated synaptic conductance reduces the effects of depolarizing θ_{NMDA} on *f* but cannot prevent the termination of activity for depolarizing enough θ_{NMDA}. Near the θ_{NMDA} value where the persistent activity ceases to exist, *f* is low even when *g*_{AMPA} is reduced (Fig. 9).

### Prediction from the model

Our analysis yields several predictions that can be tested in experimental systems such as cortical slices. First, gradual blockade of AMPA-mediated excitation is expected to decrease the number of spikes in a burst *N*_{s} and increases the bursting frequency *f*, or even transfer a synchronized bursting state into a tonic, asynchronous state. It will be interesting to look for weakly synchronized states when *g*_{AMPA} is blocked. Second, *N*_{s} is decreased when the persistent Na^{+} current is blocked. Third, gradual blockade of NMDA-mediated excitation is expected to decrease *f* and but not to change *N*_{s} considerably. Fourth, raising [Mg^{2+}]_{o} gradually will decrease *f* and increase *N*_{s}, until eventually the activity will stop. Fifth, the excitatory-to-inhibitory conductance in the somatosensory cortex includes an NMDA-mediated component.

### Functional implications

This work shows that slow NMDA synapses can support persistent synchronized bursts as long as they are not highly blocked by Mg^{2+} ions. This mechanism is relevant for epileptiform activity (Telfeian and Connors 1999). It may also be relevant for cortical tissues in vivo provided that the Mg^{2+} blockade is weak, namely θ_{NMDA} is hyperpolarized (Fleidervish et al. 1998) or the neuron is depolarized, such as during the “up state” (McCormick 2005; Steriade et al. 1993). Models of such networks with weak Mg^{2+} blockade of NMDA-mediated excitation were shown to exhibit persistent asynchronized states (Compte et al. 2000). We show that under the certain intrinsic and synaptic conditions we define, persistent synchronized bursting state states (Cardoso de Oliveira et al. 2001) can also occur in such networks.

Intrinsic bursters in cortex appear in layer V, and the low [Mg^{2+}]_{o} waves propagate in layer V alone but not in cortical slices without layer V (Flint and Connors 1996; Silva et al. 1991). The model shows that the intrinsic bursting property is important (although not necessary) for obtaining robust persistent network bursting activity at low [Mg^{2+}]_{o}. Therefore the cellular properties of layer V neurons together with the results of this model, taken together, can explain why this type of epileptiform propagates only in layer V. Other factors, such as the existence of long-range cortico-cortical connections in layer V (Connors and Amitai 1995), can also support those waves.

NMDA receptors underlying intracortical connections in young rats (P13–P21) are composed of the NR2B subunit, and the decay rate of their EPSCs is ∼108 ms (Kumar and Huguenard 2003). The decay rate of EPSCs of the NR2A subunits is faster, ∼65 ms. The NR2A subunit becomes abundant in cortex with respect to the NR2B subunit as a result of development (Quinlan et al. 1999) or rule learning (Quinlan et al. 2004). Our model predicts that NMDA-mediated persistent synchronized activity at ∼10 Hz can barely exist if most of the NMDA receptors are composed of the fast NR2A subunits. Faster synchronized activity, however, such that is found in vivo (Cardoso de Oliveira et al. 2001), can be supported even by NMDA receptors with NR2A subunits.

## APPENDIX A

### Model equations and parameters

There are *N* excitatory cells and *N* inhibitory cells, representing a spatial discretization of the continuous integro-differential equations (*Eqs. 1* and *2*). The position of the *i*th neuron is *x _{i} = i*/ρ. We specify here our reference parameter set, that is used throughout the paper unless stated otherwise.

### Excitatory cells

We use the Hodgkin-Huxley-like formulation for the excitatory cortical cells based on a previous model (Golomb and Amitai 1997).

#### Current balance equation

(A1)

#### Intrinsic Currents.

##### Sodium current I_{Na}

(A2) (A3) (A4) (A5) (A6) *g*_{Na} = 35 mS/cm^{2}, *V*_{Na} = 55 mV, θ_{m} = − 30 mV, σ_{m} = 9.5 mV, θ_{h} =−45 mV, σ_{h} = −7 mV, θ_{th} = −40.5 mV, σ_{th} = −6 mV.

##### Persistent sodium current *I*_{NaP}

(A7) (A8) .

##### Delayed rectifier potassium current *I*_{Kdr}

(A9) (A10) (A11) (A12) *g*_{Kdr} = 3 mS/cm^{2}, *V*_{K} = −90 mV, θ_{n} = −33 mV, σ_{n} = 10 mV, θ_{tn} = −27 mV,σ_{tn} = −15 mV. The kinetics *I*_{Na} and *I*_{Kdr} here are faster than those of Golomb and Amitai (1997).

##### Slow potassium current I_{K-slow}

(A13) (A14) (A15) .

##### Leak current I_{L}

(A16) .

#### Synaptic Currents.

##### AMPA current I_{AMPA}

(A17) (A18) (A19) *g*_{AMPA}=.

##### NMDA current I_{NMDA}

(A20) (A21) (A22) (A23) (Kumar and Huguenard 2003). The value of θ_{NMDA} is −*∞* for , and increases logarithmically with [Mg^{2+}]_{o} (Jahr and Stevens 1990).

### GABA_{A} current I_{GABAA}^{IE}

A24 (for disinhibited networks) or 0.05 mS/cm^{2} (for networks with inhibition), .

### Inhibitory cells

We use the Wang-Buzsáki model (Wang and Buzsáki 1996).

#### Current Balance Equation

(A25) .

#### Intrinsic Currents.

##### Sodium current I_{Na}^{I}

(A26) (A27) (A28) (A29) (A30) (A31) (A32) .

#### Delayed Rectifier Potassium Current *I*_{Kdr}^{I}

(A33) (A34) (A35) (A36) .

### Leak Current I_{L}^{I}

(A37) .

#### Synaptic Currents.

##### AMPA current I_{AMPA}^{EI}

(A38) .

##### NMDA current I_{NMDA}^{EI}

(A39) .

##### GABA_{A} variables

(A40) .

### Architecture

The excitatory-to-excitatory synaptic footprint shape (*Eqs. 10, A17*, and *A20*) is used with the discrete function *w*(*j*), where (A41) The basic unit length is the footprint length , corresponding to a typical width of a cortical column, ∼0.5 mm. Values of *L*/σ vary from 16 to 32, and ρ varies from 8 to 64. The excitatory-to-inhibitory synaptic footprint shape (*Eqs. A38* and *A39*) is . The inhibitory-to-excitatory synaptic footprint shape (*Eq. A24*) is *w*^{IE} (*j*) =] with (Golomb and Ermentrout 2002).

### Numerical methods

Simulations were performed using the fourth-order Runge-Kutta method with time step . Varying ρ above 8 has little effect on the results. Synaptic fields (*Eqs. A17, A20, A24, A38*, and *A39*) were computed using fast Fourier transform. Bifurcation diagrams (Figs. 8 and 10) were calculated with the software package XPPAUT (Ermentrout 2002).

## APPENDIX B

### Definition of local synchrony measure

The synchrony measure χ quantifies the normalized average voltage fluctuations, where the average is taken over a local population of excitatory neurons within a certain distance of order *σ* from a specific neuron. Here, we choose this neuron to be in the middle of the chain and compute averages over all the neurons within a distance *σ* from it. The population average voltage is (B1) The variance of the time fluctuation of *V*(*t*) is (B2) where *dt* . . . denotes time average over a large time, *T _{m}*. After normalization of σ

_{V}to the average over the population of the single-cell membrane potentials (B3) we define a synchrony measure χ(ρ) for the activity of a system with cell density ρ by (B4) This synchrony measure, χ(ρ), is between 0 and 1. In the limit σρ→∞ it behaves as (B5) where

*a*is a constant. In particular, if the system is fully synchronized [i.e., for all

*i*], and if the state of the system is asynchronous. Asynchronous states are unambiguously characterized only when the number σρ of neurons is infinite (Golomb and Hansel 2000). To test the level of synchrony of a neuronal system in simulations, one needs to calculate χ(ρ) for various values of ρ and determine the asymptotic value for very large σρ (Hansel and Sompolinsky 1996).

## GRANTS

The research was supported by United States-Israel Binational Science Foundation Grant 2003019 to D. Golomb. G. B. Ermentrout is funded by National Science Foundation and National Institute on Mental Health.

## Acknowledgments

We thank Y. Grossman for helpful discussions and Y. Amitai and D. Hansel for careful reading of the manuscript.

## Footnotes

↵1 As a control parameter varies, such bursting states are obtained from tonic, periodic state via a period doubling or a torus (Hopf of periodics) bifurcation (Mandelblat et al. 2001).

↵2 The auxiliary variables

*x*_{NMDA}is used just for determining*s*_{AMPA}, and is not included in the analysis.The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2006 by the American Physiological Society