INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

In the preceding paper (Butera et al. 1999), we presented minimal biophysical models of oscillatory bursting of respiratory neurons in the mammalian pre-Bötzinger complex (pre-BötC); current evidence suggests that these neurons, coupled via excitatory amino-acid (EAA)-mediated synaptic connections, form a bilaterally distributed, hybrid pacemaker network that is the origin of rhythm generation in vitro. The bursting frequency of these neurons is voltage dependent, potentially regulated by a variety of intrinsic and synaptic mechanisms that control membrane potential, including synaptic input from a tonically firing population of beating neurons (see Smith 1997; Smith et al. 1995 for review). The dynamic interactions of both the intrinsic membrane and synaptic mechanisms underlie rhythm generation. In the present paper, we investigate these dynamics in a model pacemaker network. We utilize one of the cellular models presented in the companion paper (Butera et al. 1999) to investigate how factors such as synaptic coupling strength and heterogeneity of cellular properties affect the oscillation frequency and synchrony of the population of coupled bursting pacemaker neurons. We begin by studying a pair of model neurons and analyze how excitatory synaptic coupling and excitatory tonic drive affect the mode of activity (silence, bursting, beating) and the frequency of bursting. We then consider a larger network consisting of a population of burst-capable neurons (neurons that exhibit intrinsic bursting at some level of depolarizing input) with parameter heterogeneity and study the synchronization and frequency control of this population with comparisons with experimental data. We also study emergent network rhythms that arise from populations of nonburst-capable neurons. Finally we consider an even larger case of a population of tonically firing cells providing input to a large network of coupled bursting neurons that provides synaptic drive to a population of follower cells. The bursting neurons and follower cells are heterogeneous, and the resulting spike-frequency histograms and voltage-clamp synaptic current measurements are related to experimental data. Preliminary reports of these modeling results have been presented in condensed form (Butera et al. 1997, 1998a,b).

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

All simulations were performed on Pentium-based UNIX/LINUX or SGI IRIX workstations. Simulations were coded in the C programming language. Numerical integration of simulations of pairs of cells were performed using the numerical integration package CVODE (Cohen and Hindmarsh 1996), available at http://netlib.cs.utk.edu/ode/cvode.tar.Z. Numerical simulations of large networks of cells were integrated using a fifth order Runge-Kutta-Fehlberg method with Cash-Karp parameters and an adaptive step size (Press et al. 1992). For final simulations, relative and absolute error tolerances were 10⁶ or smaller for all state variables.

Automated determination of modes of activity (silence, bursting or beating) and a quantification of these dynamics (burst duration, burst frequency, etc.) for simulations of pairs of cells required a high degree of accuracy (see DISCUSSION), thus we adopted the following procedure:

	1
	Run the simulation for 60 s (of simulation time) to allow transients to decay.
	2
	Collect interspike intervals (ISIs) for 60 s. Let p equal the maximum of these ISIs. Let w define a window size equal to 0.9p.
	3
	If there were no spikes, the cell is classified as silent.
	4
	Collect ISIs for 120 s. Collect the following statistics: ISI, burst duration (BD or active phase duration), burst interval (BI or silent phase duration), and burst period (BP). A sliding window of size w is used to determine the transition between active and silent phases of a burst cycle. The beginning of a burst is defined at any spike following an interval of duration w or longer. The end of a burst is defined at the last spike preceding an interval of duration w or longer.
	5
	If over the collection period there was never a transition to a silent phase, the cell is classified as beating and the mean, standard deviation, minimum, and maximum of the ISIs is computed.
	6
	Otherwise, record the mean, standard deviation, minimum, and maximum of the ISIs, BDs, BIs, and BPs. Let c equal the coefficient of variation (standard deviation divided by the mean) for both the ISIs and BPs. If cISI < s × cBP, the cell is classified as beating. Otherwise, the cell is classified as bursting. The value of s was determined empirically, and for the present results, we used s = 0.2. We chose a value of s < 1 because in the case of coupled populations beating solutions tend to show little variability, while solutions that we would visually classify as bursting still may have an inherent degree of variability in burst period from cycle to cycle. Qualitatively similar results were obtained for values of s from 0.1 to 0.95.

The concept behind the above algorithm is that bursts may be identified from spike trains by looking for particularly large ISIs. Thus ISIs that are "similarly large" (the window size 0.9w) are used to identify silent phases between bursts. For simulations of large populations (50-500) of bursting neurons, we adopted a simpler definition of the transition between the active and silent phase of the burst. A burst begins at the first spike following a window of 500 ms. From these spike times, the BD, BI, and BP are calculated. For long time series of spike times, these values were calculated for each cycle and averaged.

Model of synaptic dynamics

Throughout this paper we will use model 1, as presented in the preceding paper (Butera et al. 1999). During oscillatory bursting activity, the burst is initiated by a persistent Na⁺ current (I_NaP-h) and terminated by slow voltage-dependent inactivation of that current. We have focused on this model because its dynamics are more consistent with our experimental recordings; however, some of the general results presented in this paper (such as those presented in Figs. 1 and 6) have been reproduced qualitatively using model 2 as well.

View larger version (28K): [in this window] [in a new window]   Fig. 1. Effect of level of depolarization (via gtonic-e) and excitatory coupling on synchronous activity of 2 coupled bursting cells. Only activity of 1 cell is shown because in all cases of bursting the bursts are synchronous. Left: nonbursting quiescent modes. Panels with neuron spiking appearing as a black band are beating states. Inset: beating activity at an expanded time scale. For gtonic-e = 0.2 and syn-e = 15 nS the coupled system is bistable, with coexisting silent and beating solutions (top left). Parameter values of gtonic-e and syn-e were chosen so as to represent the range of dynamic behaviors that are exhibited by the coupled system.

The equations and parameter sets are those presented in the companion paper (Butera et al. 1999) unless otherwise noted. For simulations of large populations, some parameters were picked randomly from a normal distribution with a given mean and standard deviation (see Table 1). To model intercellular coupling, we have added an additional current, I_syn-e, to the model. Thus the equation for the rate of change of the membrane potential is modified to be

(1)

The pacemaker neurons in the pre-BötC are postulated to receive both excitatory and inhibitory input from populations of tonically firing (beating) cells (Funk et al. 1993; Greer et al. 1991; Smith et al. 1995). For simplicity we only consider a mean steady level of excitatory input, represented by the conductance g_tonic-e, although analogous results are obtained if a mean level of inhibitory input is considered as well. Synaptic input from this beating population is modeled as I_tonic-e = g_tonic-e(V  E_syn-e), where E_syn-e is 0 mV, the reversal potential of non-N-methyl-D-aspartate (non-NMDA) EAA-mediated synaptic currents, the type of EAA current thought to mediate tonic excitation (Funk et al. 1993; Smith et al. 1991). This current is used in all simulations (Figs. 1-3) where only the mean level of postsynaptic conductance is included and we do not explicitly model the population of tonically beating neurons, otherwise g_tonic-e = 0 nS.

I_syn-e models the EAA-mediated coupling between individual bursting neurons in all simulations. It also is used in those simulations that explicitly consider a modeled population of beating neurons, in lieu of I_tonic-e. The synaptic input to neuron j from the population of N neurons is described as

(2)

where the gating variable s is driven by presynaptic depolarization according to

(3)

This formulation is based on the assumption that the binding of transmitter is rapid compared with the rate of the conformational change of the channel (Magleby and Stevens 1972) and is similar to the formulation described in Wang and Rinzel (1992). The postsynaptic responses (_s = 5 ms, k_r = 1) are assumed to be fast due to non-NMDA glutamatergic receptor activation. The function s(V_pre) is a sigmoidal function whose half-activation value is set such that it is only activated by the firing of a presynaptic action potential. The function s(V) is 1/{1 + exp[(V  _s)/_s]}, where _s = 10 mV and _s = 5 mV. For population simulations, _syn-e was chosen randomly for each cell from a normal distribution with a specified mean and standard deviation (Table 2).

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Pairs of cells

To study the effects of excitatory coupling on the dynamics of bursting, we initially consider the simple case of a pair of neurons. Each neuron of the pair uses the ionic current equations and parameters of model 1 (Butera et al. 1999). These neurons are coupled reciprocally by an excitatory synaptic current (I_syn-e). Each neuron in the pair also receives a mean level of tonic excitatory synaptic drive (I_tonic-e). As described in the companion paper, this current has a similar depolarizing effect on the burst dynamics of single neurons as I_app. We studied the dynamics of this pair of bursters as _syn-e and g_tonic-e were varied to assess the control of bursting frequency by excitatory tonic drive (g_tonic-e) and excitatory synaptic coupling (_syn-e). In all cases, g_tonic-e and _syn-e were identical for each cell. The time course of the membrane potential oscillations for various values of g_tonic-e and _syn-e are illustrated in Fig. 1. The frequencies of bursting and beating, as well as burst duration and burst period, at various values of _syn-e as g_tonic-e is varied are plotted in Fig. 2.

View larger version (44K): [in this window] [in a new window]   Fig. 2. Effects of gtonic-e and syn-e on the dynamics of bursting and beating of pairs of coupled cells. Plots are of burst duration (A), burst interval (B), burst period (C), duty cycle (D), beat frequency (E), and burst frequency (F). Colored areas indicate regions of dynamic activity being quantified. Range of shaded values for each diagram is indicated by the associated range bar for each figure. Units are indicated at the top of each range bar. Color values were scaled to parameter ranges that most accurately revealed trends in the data, and values beyond the minimum or maximum of the range of values used for coloring were clipped at those values.

Both g_tonic-e and _syn-e activate excitatory currents, although g_tonic-e may be described as constant and _syn-e as phasic. The cases of weak coupling (_syn-e < 1 nS) and high tonic drive (g_tonic-e > 0.6 nS) will be considered shortly. For the rest of the parameter space illustrated in Fig. 2 where bursting activity occurs, several trends are evident.

At a given level of g_tonic-e, _syn-e increases both the burst duration (Fig. 2A) and the burst interval (Fig. 2B), with a net effect of increasing the burst period. These effects are not independent: increasing excitatory coupling, which is triggered by the firing of action potentials, causes action potentials in each cell to further excite the other cell, with a net increase in intensity of firing in both cells. The higher spike frequency increases excitatory synaptic input during the active phase of a burst. For the burst to terminate, I_NaP-h must be inactivated further than it would in the absence of synaptic input. Thus the burst duration is extended. The additional inactivation of I_NaP-h prolongs recovery from inactivation as well, increasing the duration of the silent phase of the burst cycle. Hence we obtain the net, counterintuitive, effect of reduced burst frequency with stronger excitatory coupling (also noted by Somers and Kopell 1993 for idealized slow-wave neural oscillators).

At a given level of _syn-e, burst duration increases with g_tonic-e. Through most of this range, the burst interval decreases with g_tonic-e. This decrease in burst interval is the dominating effect and the overall burst period decreases. However, at higher values of g_tonic-e in this parameter range, the burst interval changes minimally or slightly increases. In this region, which is just before the region of complex dynamics described later in this paper, the increase in burst duration is dominant and the burst period increases.

An important effect of synaptic coupling (except possibly for very strong _syn-e) is that it extends the range of values of g_tonic-e over which bursting is supported. This effect will be investigated further in the following two sections, which consider the dynamics of larger populations of bursting neurons. This increased dynamic range is maximal at approximately _syn-e = 2 nS (Fig. 2). The half-pear shaped regions extend rightward beyond their intersection with the abscissa near g_tonic-e = 0.4 nS, which denotes the bursting-to-beating transition point for an isolated cell. The transition from silence to bursting is the vertical boundary (on left) because there is a critical value of g_tonic-e necessary to depolarize the cell into a bursting mode, and _syn-e is inactive when both cells are silent. In contrast, the transition between bursting and beating occurs at a voltage range where the coupling synapses are fully turned on, and there is an approximate trade-off between g_tonic-e and _syn-e that determines the transition, i.e., a critical amount of synaptic input (both tonic or excitatory coupling) is necessary to maintain the cells in a beating mode. Less g_tonic-e requires more _syn-e to maintain beating and vice versa. This explains the nearly linear slope of the bursting regime's right-side boundary, above the regime of weak coupling. This trend is only applicable at high coupling strengths (>1 nS).

At low coupling strengths (_syn-e < 1 nS), synaptic input during the firing of action potentials synchronizes the two cells but does not significantly alter the dynamics of each cell (compared with the uncoupled state). In this case, the bursts of the two cells synchronize but the burst duration hardly changes with g_tonic-e, while the burst interval decreases, as in the case of the uncoupled cell (_syn-e = 0); see Fig. 2, A-C. The net effect on burst period is similar to that in the case of stronger coupling: an overall decrease in the burst period. Also, unlike the case of strong coupling strengths (_syn-e > 1), increasing _syn-e increases the value of g_tonic-e where the transition from bursting and beating takes place.

At high values of g_tonic-e (>0.6), regardless of the value of _syn-e, the model neurons approach the parameter range at the interface between bursting and beating and their dynamics become quite complicated and sometimes irregular. This is most evident in the burst duration (Fig. 2A), and those quantities that depend on this measure (Fig. 2, C and F).

Effects of coupling on population burst dynamics

Real neurons (as opposed to deterministic model neurons) possess considerable variability in their intrinsic membrane parameters, such as maximal conductances. We investigated the role of coupling in synchronizing a population of bursting neurons with heterogeneous properties. We chose as our heterogeneity parameters E_L, _NaP-h, and _syn-e. E_L determines the intrinsic baseline level of depolarization of the bursting cells, an important parameter for voltage-dependent frequency control as described in Butera et al. (1999); _NaP-h determines the ability of a neuron to intrinsically burst (see Figs. 6 and 7, Butera et al. 1999); and _syn-e is the coupling strength between individual neurons. Tables 1 and 2 specify the intrinsic and synaptic parameters that were allowed to vary and how they were distributed for all of the population simulation results presented. Typically parameters were chosen randomly for each cell from a normal distribution with a specified mean and standard deviation. The standard deviations of E_L and _NaP-h were chosen so that the distributed parameters fell within the range of dynamic behaviors studied in Butera et al. (1999). For the random distribution of synaptic conductances, we specified a standard deviation of 25% (or more) of the mean value. A condition was enforced that all conductances be greater than zero. For all simulations, initial conditions were chosen randomly from a uniform distribution across the physiological dynamic range of individual state variables. After any change in parameters or initial values, a settling period of 60 s of simulation time was allowed before data were collected.

Some of the simulations presented in this section were repeated where all the intrinsic conductances of the model were selected randomly from a normal distribution, using the nominal conductances for model 1 in Butera et al. (1999) as the mean and a CV of 20%. Results were similar to those presented in this section, where the intrinsic heterogeneity parameters are solely E_L and _NaP-h.

Figure 3 illustrates the dynamics of individual neurons and overall network activity before and after synaptic coupling for a typical simulation. The population consisted of 50 heterogeneous bursting neurons with all-to-all coupling. The left column illustrates the dynamics of the uncoupled network (_syn-e = 0), and the right column illustrates the dynamics of the coupled network. When uncoupled, neurons in the heterogeneous population exhibit a variety of activity modes, with some spiking tonically, some bursting, and some remaining silent (raster plot, Fig. 3A1). There is little correlated network activity (population spike histogram, Fig. 3A2). This dispersion in intrinsic activities is consistent across 50 simulations with different randomly generated parameter sets and initial conditions (bar plots, Fig. 3A3). After coupling, most of the cells burst in relative synchrony with a coordinated periodic burst of action potentials (Fig. 3B, 1 and 2). This finding is also consistent across a large number of simulations (B3).

View larger version (41K): [in this window] [in a new window]   Fig. 3. Effects of excitatory amino acid (EAA) coupling on synchronization of a heterogeneous population of 50 bursting neurons. Raster plots illustrate spike times of the population of the 50 cells before (A1) and after (B1) coupling. Network activity (A2 and B2) is defined by calculating histograms of spike times across the population (bin size = 10 ms). Bar plots illustrate number of cells (means ± SD) in population whose dynamics are defined as silent, bursting, or beating before (A3) and after (B3) coupling. Results are pooled from 50 simulations.

Among the population of burst-capable neurons, the excitatory coupling allows the bursting cells to recruit silent neurons by providing sufficient additional depolarization to trigger some of the silent cells to burst. Tonically spiking cells also are recruited into bursting via additional excitatory input from bursting cells, which transiently increases spike frequency. This increase in frequency leads to additional inactivation of I_NaP-h, which temporarily may terminate spiking, resulting in bursting behavior. The recruitment of both intrinsically silent and spiking cells to a bursting mode increases the net excitatory synaptic activity in the population, making bursting at a network level more robust (see following text).

The relatively synchronous burst dynamics were found to be quite tolerant of parameter heterogeneity. We performed numerous additional simulations and report several initial observations regarding the effects of coupling and parameter heterogeneity on the dynamics of the network. First, general synchrony of bursts is obtained if the coupling is sufficiently strong. Greater parameter heterogeneity requires stronger coupling to synchronize bursts (Pinsky 1994; Pinsky and Rinzel 1994). Second, we investigated the effects of sparse coupling, where each cell-to-cell connection has a probability of existing. We found that connectivity patterns did not greatly affect the ability of the network to synchronize bursts (even with connection probabilities as low as 3%), given a similar level of mean synaptic input to each cell. However, although sparsity did not greatly affect synchronization of the population, the interval between burst episodes of the population became more irregular as the connectivity was made sparser. Third, when synchrony starts to break down, due to increased parameter heterogeneity, weaker coupling, or by decreasing the mean level of E_L such that less cells are intrinsically bursting or spiking, the network usually makes a transition to a mode where most cells burst every N cycles, with some bursting at higher multiple frequencies (Fig. 4, A and B). In this case, the N:1 phase-locking of the cell to the network output is still quite regular. In some cases, hyperpolarizing the mean level of E_L resulted in irregular network-wide synchronized bursts (Fig. 4C).

View larger version (35K): [in this window] [in a new window]   Fig. 4. Network activity modes with partial synchrony. A1-C1: raster plots. A2-C2: network activity, defined as histograms of spike times across the population (bin size = 10 ms). A, 1-2, and B, 1-2: 2:1 and 4:1 synchrony, where a fully synchronized burst across the entire population of cells occurs every 2 or 4 cycles, respectively, with some cells bursting at more frequent intervals, resulting in partial population synchrony with lower amplitude population bursts. C: irregular burst firing with varying degrees of synchrony resulting in both high- and low-amplitude population bursts.

Johnson et al. (1994) performed a systematic study of the effects of synaptic coupling on the dynamics of oscillatory bursting neurons in the pre-BötC. Inspiratory neurons in the pre-BötC were identified via extracellular recording, and those that continued to burst rhythmically after switching to a low-Ca²⁺ solution to block synaptic transmission were classified as possessing intrinsic bursting properties. The burst duration and burst frequency were quantified for each bursting cell in the control solution (coupled network) and low-Ca²⁺ solution (uncoupled), allowing an assessment of the role of coupling on the dynamics of bursting at a cellular level. Figure 5 is a replotted version of the data presented in Fig. 8, A and B, from Johnson et al. (1994). The general effects of excitatory coupling on burst frequency are illustrated in the histogram: the cell-to-cell variance in burst frequency is reduced, and the distribution of burst frequencies is compressed toward lower (slower) values. The mean burst duration increases with coupling. These effects are similar to those illustrated in the simulations with only two cells (Fig. 1).

View larger version (31K): [in this window] [in a new window]   Fig. 5. Comparison of burst parameters for in vitro bursting inspiratory neurons under low-Ca2+ conditions () and in control solution where neuron bursting was synchronous with respiratory network rhythm (). Scatter plots illustrate burst frequency (A) and burst duration (B) for each neuron in a control and low-Ca2+ solution. , line of identity. Histograms summarize total number of neurons with parameter values falling within bins indicated. Data taken and replotted from Fig. 8 of Johnson et al. (1994) with permission.

We performed numerous independent simulations of our network of 50 bursting cells using the heterogeneity parameters specified in Table 3. These simulations were run without and with EAA coupling as previously illustrated in Fig. 3. For each cell in each simulation that burst when uncoupled, we quantified the burst frequency and burst duration before and after EAA coupling. This yielded a large number of cells on which to quantify the effects of coupling in a heterogeneous network. The pooled results (Tables 4 and 5), when quantified, show similar effects as those in Fig. 5. However, Johnson et al. (1994) typically only recorded from one cell per slice preparation, thus having quantifiable data from <50 bursting neurons. To assess the role of sampling in potentially biasing the results, we randomly selected only one endogenously bursting neuron from each of 50 simulations and pooled the results, which are illustrated in Fig. 6. Similar trends are evident even in the much reduced dataset: coupling decreases the mean and standard deviation of the burst frequency, while coupling increases the mean burst duration (see Table 5).

View larger version (29K): [in this window] [in a new window]   Fig. 6. Comparison of burst parameters for in computo bursting inspiratory neurons without () and with EAA coupling where neuron bursting was synchronous with respiratory network rhythm (). Data were compiled by selecting 1 bursting neuron at random from each of 50 simulations to approximate neuron sampling using the protocol specified in Johnson et al. (1994). Scatter plots illustrate burst frequency (A1) and burst duration (B1) for each neuron in a coupled and synaptically uncoupled population. , line of identity. Histograms (A2 and B2) summarize total number of neurons with parameter values falling within bins indicated.

We used various other parameter distributions (Table 3) for additional simulations to assess the robustness of these coupling effects. We chose different parameter distributions to vary the number of intrinsically silent, bursting, and beating neurons (Table 4). These distributions ranged from a nearly even distribution of silent, bursting, and beating neurons (set 1) to cases with a majority of silent neurons and <20% intrinsically bursting (sets 6 and 7). In general, we found that the coupling effects described earlier persisted across all the parameter distributions where synchronous bursting activity occurred (Table 5). In Table 5, set 1 corresponds to the results presented in Figs. 3, and set 1r corresponds to the results presented in Fig. 6.

The results obtained in this section were verified further by running simulations of a population of 500 cells using a similar parameter distribution as in parameter set 1 with 1/10 of the cell-to-cell coupling strength (_syn-e) because there are 10 times as many cells. The parameters for this simulation are referred to as parameter set 8 (Table 3), and the effects of coupling on activity mode and burst dynamics are specified in Tables 4 and 5, respectively. Results are similar to those obtained for parameter set 1, thus the results reported in this section using a population of 50 cells do not appear to be altered by small population effects.

Frequency control of network activity

As analyzed in the companion paper (Butera et al. 1999), the oscillatory bursting neurons in the pre-BötC exhibit voltage-dependent frequency control, where the burst frequency of individual neurons is regulated by depolarizing input, whether intrinsic (E_L) or synaptic (g_tonic-e). We investigated how well this mechanism of frequency control persisted in a large population of neurons with the same types of parameter heterogeneity described in the previous section. For each simulation, the mean level of E_L was set to 65 mV. The mean depolarizing input to the population, g_tonic-e was initially zero and increased every 120 s of simulation time in 0.05-nS increments. Figure 7A, 1-7, illustrates the aggregate population spike activity of one simulation as a function of g_tonic-e. Parameters were distributed as indicated in Tables 1 and 2.

View larger version (34K): [in this window] [in a new window]   Fig. 7. Voltage-dependent frequency control of network activity. Left: histograms of network spiking (10-ms bin size) as the mean level of gtonic-e is increased. B: network burst frequency as a function of gtonic-e. C: burst frequency of a single bursting model cell using the mean parameter values of the network simulation.

As g_tonic-e is increased (Fig. 7A, 1-6), the following effects on population bursting are exhibited: an increase in the population burst frequency; a decrease in the amplitude (spike frequency) of the burst of network activity; an increase in the spread of the onset of burst firing (the bursts of network activity at low g_tonic-e have sharp rises/falls and large peak amplitudes, while the bursts of network activity at higher levels of g_tonic-e have more gradual rises/falls and smaller peak amplitudes); and a decrease in the signal to noise ratio, as evidenced by an increase in the mean and variance of spike activity and a decrease in the amplitude of the network burst. At a sufficient level of depolarization (Fig. 7A7), all of the neurons in the network are spiking continuously with no coordinated bursting evident. A comparison of the burst frequency of the network simulation and the range of burst frequencies of a single cell using the mean values of E_L and _NaP-h and identical g_tonic-e values is illustrated in Fig. 7C. This reveals that the dynamic range of g_tonic-e where bursting is supported is much larger for the population than for a single cell. Similar effects are obtained in the two cell case (Fig. 2).

To assess which factors contributed to the g_tonic-e versus burst frequency relationship of the heterogeneous network, we ran additional network simulations with no parameter heterogeneity (all cells have identical parameters), intrinsic heterogeneity only (E_L and _NaP-h), and synaptic heterogeneity only. The effects of increasing coupling strength also were investigated. The results are illustrated in Fig. 8. The homogeneous network had a dynamic range of g_tonic-e larger than the single cell but less than that for the heterogeneous network (Fig. 8A). However, the homogeneous network displayed a smaller range of burst frequencies than both the single cell and the heterogeneous population. We also found that synaptic parameter heterogeneity yielded a frequency versus g_tonic-e relationship similar to the homogeneous network, whereas intrinsic parameter heterogeneity yielded a frequency versus g_tonic-e relationship similar to the fully heterogeneous network (Fig. 8B). From these results, we conclude that both intrinsic heterogeneity and excitatory synaptic connectivity contribute toward the increased dynamic range of g_tonic-e where bursting occurs (as opposed to a single mean-value cell); intrinsic heterogeneity is necessary for the population to have a range of burst frequencies similar to that of the average single cell response; and synaptic heterogeneity made little difference as long as the coupling did not become effectively too sparse.

View larger version (19K): [in this window] [in a new window]   Fig. 8. Effects of parameter heterogeneity and coupling strength on frequency control of network bursting. Curves plot burst frequency as a function of gtonic-e. A: heterogeneity effects. Comparison of ranges of burst frequencies for the single-cell, homogeneous network, and heterogeneous network. The heterogeneous network and single-cell curves are identical to Fig. 7, B and C, respectively. B: sources of heterogeneity. Comparison of frequency control of the heterogeneous network (EL, NaP-h, syn-e) with networks with intrinsic heterogeneity only (EL, NaP-h) and synaptic heterogeneity only (syn-e). Both of these networks used the exact same parameter distributions as the heterogeneous network for the appropriate heterogeneity parameters. Curve for intrinsic heterogeneity is superimposed partially on curve for heterogeneous network. C: coupling effects. Comparison of frequency control of the heterogeneous network with identical networks where the coupling strengths among the bursting cells have been doubled (×2) or tripled (×3).

We also investigated the role of coupling strength on the burst dynamics of the network (Fig. 8C). The primary effect is similar to that observed in the two-cell simulation (Fig. 1): increasing the coupling strength increases the mean burst duration and decreases the mean burst frequency of the population. As the coupling strength was increased, the ranges of both g_tonic-e and burst frequencies where synchronous bursting occurred were decreased. We also found that higher coupling strengths also decrease the temporal distribution of burst onset times in heterogeneous populations (not shown). At very high coupling strengths, the network showed only beating activity with no subthreshold oscillations at the cellular level, analogous to the phenomena known as "oscillator death" (Ermentrout and Kopell 1990). This differs from the case of modeled bursting neurons with gap-junctional coupling (Sherman and Rinzel 1992), where high coupling strengths lead to the dynamics of the coupled cells approaching that of a "mean" cell.

We repeated the simulations of this section using the mean value of E_L, instead of g_tonic-e, as a frequency-control parameter. Similar results to those shown in Figs. 7 and 8 were obtained.

As g_tonic-e is varied through the range where the network bursts rhythmically, the distribution of intrinsic cell firing properties changes (Fig. 9). At the low end, coupling is adequate to induce bursting even though nearly 90% of the cells are silent. At the upper end, bursting is maintained with a population that is >90% beaters. This graphic illustration emphasizes that this population is made up of voltage-dependent (conditional) bursters: cells with I_NaP-h that can be recruited readily to burst by collective depolarizing inputs even though few of the cells are spontaneous bursters when decoupled from each other. For each g_tonic-e within the network's operational range, the synchronized population burst frequency is below the mean frequency of the individuals. Again, this behavior is the large-population analogue (including heterogeneity) for the counterintuitive effect of phasic excitatory synaptic input on burst frequency noted for a pair of identical cells in Fig. 2. Figure 9 also can be used to predict the results of a synaptic blocking experiment in which both interneuronal coupling among the bursting population as well as excitatory tonic drive are removed (see DISCUSSION).

View larger version (35K): [in this window] [in a new window]   Fig. 9. Intrinsic modes of activity for disconnected component cells in simulation of Fig. 7. A: each horizontal line represents the range of gtonic-e where individual cells exhibit intrinsic bursting. B: number of cells demonstrating intrinsic silence, bursting, or beating activity modes as a function of gtonic-e. C: , mean, minimum, and maximum burst frequency of cells with intrinsic bursting activity as a function of gtonic-e. For comparison, , synchronous burst frequency as a function of gtonic-e for the coupled network.

Emergent network oscillations from coupled cells with I_NaP-h

The examples above have used populations of burst-capable neurons where at any given parameter set, some fraction of the neurons were intrinsically bursting. Simulations presented in this section examined network behavior with a population of 50 neurons with a low mean level of _NaP-h such that no neuron in the population exhibited intrinsic oscillatory bursting properties at any level of depolarization (E_L or I_app), i.e., they were only capable of spiking or silence. In the companion paper (see Fig. 7 in Butera et al. 1999) we define the parameter space in _NaP-h and E_L determining the activity modes of the pacemaker cells, which indicates regions of _NaP-h with only quiescent or beating modes. We therefore wanted to determine if there were low _NaP-h parameter sets where synchronous oscillation could emerge as a network property. In the simulation shown in Fig. 10A, the uncoupled population shows only five neurons exhibiting spiking activity, the rest are silent. Coupling the population (Fig. 10B) increases the level of spiking activity and recruited two more neurons to spiking, but the rest of the population remains silent. An additional 25% increase in _syn-e gave rise to synchronous bursting (Fig. 10C) across the entire population. Alternatively, an increase in the mean value of _NaP-h (Fig. 10E) also gave rise to synchronous population bursting. Even at this higher mean value of _NaP-h, none of the neurons in the population exhibit intrinsic bursting activity (Fig. 10D). In these cases (Fig. 10, B and E), a sufficient amount of depolarizing subthreshold inward current from both intrinsic (_NaP-h) and synaptic (_syn-e) sources will initiate a burst, and a decrease in one current can be offset by an increase in the other. We identified two other necessary criteria for synchronous bursting activity to occur when _NaP-h is low: 10% or so of the population must possess intrinsic spiking activity and sufficient _NaP-h must be present to inactivate and terminate the active phase of the burst. In effect, the combined actions of _syn-e (synaptic) and _NaP-h (intrinsic) depolarize the cells to a spiking state, and the spiking recruits inactivation of I_NaP-h resulting in a transient pause in the spiking. The role of slow processes (e.g., adaptation) in generating network bursting behavior has been investigated in the case of networks of spiking neurons with adaptation currents and excitatory coupling (van Vreeswijk and Hansel 1998).

View larger version (19K): [in this window] [in a new window]   Fig. 10. Emergent synchronized bursting from a population of neurons that do not possess intrinsic bursting properties. Mean values of NaP-h and syn-e indicated on raster plots for each simulation of 50 cells: standard deviation of distributions is indicated in Tables 1 and 2. For these particular distributions of intrinsic parameters (A and D), none of the neurons in the population exhibited endogeneous bursting behavior at any value of EL or Iapp. Synchronized bursting and network oscillations emerge (C and E) with temporal dispersion in the onset of spiking of different neurons.

The oscillations that occur in these cases represent an emergent property of the network of excitatory coupled neurons with low _NaP-h. The oscillatory frequencies generated under these conditions tend to be slower and at the lower end of those achieved where _NaP-h is high enough for intrinsic burst generation. Preliminary simulations of frequency control of the network under low _NaP-h via g_tonic-e (not shown) reveal a restricted range of frequencies where bursting occurs, similar to the case of strong coupling in Fig. 8C (×3).

Components of network activity

In this section, we investigate the dynamics of a more complete model network that incorporates the pacemaker population kernel, a population of beating neurons that regulate the excitability of the kernel, and a follower cell population that is synaptically driven by the pacemaker cells. These populations are hypothesized to be the rudimentary elements generating, controlling, and transmitting the rhythm in the pre-BötC in the in vitro slice preparation (Smith et al. 1991, 1995). The model network consists of the following components:

	1
	50 bursting neurons with all-to-all EAA-mediated coupling. The parameters NaP-h, EL, and syn-e are randomly distributed (Table 1), similar to the simulations of the preceding section.
	2
	200 beating neurons, which represent the population of beating cells providing tonic input to the bursting neurons (in lieu of the mean tonic input parameter gtonic-e used in earlier simulations). For simplicity, we use our bursting cell model with NaP-h set to zero, resulting in a model that is capable only of silent or beating modes of activity. The frequency of spiking is determined by EL. EL was distributed among the beating cells as indicated in Table 1. This distribution of EL yields a population of beating neurons with 97% of the frequencies in the range of 1-15 Hz with a mean of ~10 Hz, consistent with spiking frequencies of tonically active neurons found in vitro (Johnson et al. 1994; Smith et al. 1990). The exact parameters of these cells are not important because they do not receive any input from within or outside the population and only serve to generate asynchronous tonic input to the population of bursting neurons. The distribution of postsynaptic conductances is indicated in Table 2. Each cell in the beating population has a 20% chance of connecting to each of the cells in the bursting population. This was done to ensure that the tonic input to individual bursting neurons was not highly correlated.
	3
	50 follower neurons that receive non-NMDA EAA excitatory input from the population of bursting neurons. The population activity of these follower neurons transmit the rhythm to (pre)motoneurons and may be considered to be indicative of integrated recordings of motor output from ventral roots of the hypoglossal nerve in the in vitro slice preparation (see Funk et al. 1993; Smith et al. 1991). We again used our bursting neuron model with NaP-h set to 0. Thus these neurons are incapable of intrinsically bursting and function as follower cells that fire action potentials when the synaptic input exceeds a certain level. The connectivity and random distributions of EL and syn-e are specified in Tables 1 and 2. Values of syn-e were chosen so that the amplitudes of the total synaptic current envelope under simulated voltage clamp were within the range observed in vitro for follower-type neurons (e.g., Funk et al. 1993; Smith et al. 1992). Each cell in the follower population has a 10% chance of receiving input from each cell in the bursting population. Thus the follower population is heterogeneous in intrinsic properties, synaptic conductances, and synaptic connectivity.

Figure 11 illustrates how these populations combine to produce rhythmic network activity. In Fig. 11A, the network is disconnected. The beating neurons are not providing input to the bursting neurons and the bursting neurons are uncoupled, some of which are firing bursts independently in an unsynchronized fashion. The raster plot shows that with the particular randomized parameter distribution used, 15 of the cells are bursting, 11 of the cells are beating, and 24 of the cells remain silent. In Fig. 11B, the bursting neurons still are uncoupled but receive depolarizing synaptic input from the population of beating cells. The bursting neurons are more excitable than in Fig. 11A and burst with less regularity due to the summated asynchronous synaptic input from the beating cell population. In this simulation, 30 neurons are bursting, 18 are beating, and 2 neurons remain silent. To assess how the effects of a "noisy" synaptic input contributed to the recruitment of bursting neurons, we reran the simulation using an identical assignment of intrinsic and synaptic parameters, only using a mean synaptic conductance (0.19 nS) via g_tonic-e, in lieu of a population of beating neurons. In this case, only 18 of the neurons were classified as bursting, with 23 beating and 9 silent. We speculate that asynchronous "noise-like" ongoing input to the population of burst-capable neurons is an additional factor that promotes bursting behavior in individual neurons. The noise occasionally will kick cells out of either the silent or beating mode to execute transient bursts. However, such a role for asynchronous synaptic input is predicated on the assumption that the number of synaptic inputs to each cell is "finite" (i.e., not a very large number).

View larger version (42K): [in this window] [in a new window]   Fig. 11. Effects of EAA synaptic input on a network of bursting neurons. Network consists of 200 beating neurons, 50 bursting neurons, and 50 follower neurons. In each panel, the 3 voltage traces show 10 s of activity in 3 bursting cells, along with the summated spike-activity (10 ms/bin) of the follower cell population (F). Raster plots (right) illustrate the spike-firing times for all 50 bursting cells for 10 s. A: bursting neurons are uncoupled from each other and the tonic cells. B: bursting neurons receive excitatory input from the tonic cells but still are uncoupled from each other. C: bursting neurons receive excitatory input both from the tonic cells and from the synaptic coupling within the bursting cells population. Vertical dashed lines illustrate reference times used for calculation of spike-frequency histograms in Figs. 12 and 13. Each reference time corresponds to a positive crossing of 20 spikes/bin of the summated spike activity of the following population. C also illustrates a raster plot (bottom) of the spike activity of the follower population (R) that is summated to provide the activity measure of the follower population (F). See text for details.

In Fig. 11C, the bursting neurons are coupled to each other and the entire population bursts in a coordinated fashion (see raster plot). The summated effect of this coordination is represented in the regular burst firing of the follower neuron population as well. The more depolarized bursting neurons (high E_L) spike earlier in the cycle, whereas the more hyperpolarized neurons burst later in the cycle. Raster plots of the spike activity of the population of follower neurons and the integrated population activity are illustrated in Fig. 11C. The main bursts of the coupled pacemaker and follower cell population activity exhibit a rapidly peaking and slowly decrementing time course characteristic of respiratory network activity in in vitro preparations (e.g., Funk et al. 1993; Smith et al. 1990).

Heterogeneity and spike-frequency histograms

Respiratory neurons conventionally have been described by their firing patterns, which are referenced to the epoch of the respiratory cycle (inspiratory, expiratory) during which the neuron fires. Johnson et al. (1994) classified the firing patterns of respiratory neurons with intrinsic bursting properties in the pre-BötC with respect to the onset of XII motor output, which represents inspiratory phase activity (Smith et al. 1990, 1991). Spike-frequency histograms were computed for the neurons using the onset of hypoglossal discharge as a reference when synaptic transmission was intact; intrinsic bursting properties of these cells were identified after blocking synaptic transmission with low-Ca²⁺ medium. A majority of the neurons (34/67) were I cells, i.e., they burst with an onset of spiking in synchrony with the XII output. A minority of cells (11/67) were classified as pre-I cells. These pre-I cells fire action potentials at a low frequency before the motor output and often increase in firing frequency immediately before the motor output event. One cell was a late-I cell, where the onset and peak firing was after the initiation XII motor output. Representative histograms of these different three types are illustrated in Fig. 12 This temporal dispersion of spiking onset has been hypothesized to result from heterogeneity of pacemaker cell intrinsic and synaptic properties (Smith et al. 1995). We analyzed spiking patterns within our heterogeneous network model.

View larger version (19K): [in this window] [in a new window]   Fig. 12. Spike-frequency histograms of intrinsically bursting respiratory pre-BötC neurons. Left: integrated XII motor output and spike-frequency histograms of pre-I, I, and late-I pre-BötC bursting neurons from Johnson et al. (1994), reproduced by permission. Right: cycle-triggered average of the follower population spiking activity (F) and the spike-frequency histograms of cells 1-3 in Fig. 11. Binwidth for histograms is 10 ms.

Figure 12 illustrates the spike-frequency histograms of the three labeled bursting neurons (1-3) and the integrated output of the follower population (F) in Fig. 11. These histograms were calculated as described in METHODS, using the activity of the follower population to generate the reference times for each burst. The shape and timing of the spike-frequency histograms are qualitatively similar to those of the data. The more depolarized cells fire before the onset of the motor output (pre-I) and increase in firing frequency immediately preceding the motor event. The raster plot in Fig. 11C indicates a substantial dispersion in the onset of low-frequency spiking of different pre-I cells with a few cells spiking throughout the interburst interval, analogous to the phase-spanning I cells (not shown in Figs. 12 and 13) described by Johnson et al. (1994). The least depolarized cells within the bursting neuron population fire in synchrony with the motor output, whereas the most hyperpolarized cells fire after a delay.

View larger version (25K): [in this window] [in a new window]   Fig. 13. Spike-frequency histograms (A, 1-3, and D, 1-3), sample membrane potential trajectory (B, 1-3, and E, 1-3), and total synaptic current (C, 1-3, and F, 1-3) under simulated voltage clamp at 60 mV from the burster (A-C) and follower (D-F) populations. Cells from the burster population (A, 1-3) correspond to those chosen for the histograms in Fig. 12B and were calculated with respect to reference times (vertical dashed line) obtained from the follower cell population activity. Voltage traces are typical for 1 burst cycle. Voltage-clamp synaptic current traces show 1 typical cycle (gray) and the synaptic current averaged across 10 burst cycles using the same reference times as the histograms. Baseline synaptic current is due to input from beating cells. Horizontal dotted lines indicate 60 mV (B, 1-3, and E, 1-3) or 0 pA (C, 1-3, and F, 1-3).

Figure 13, A-C, illustrates typical membrane potential trajectories (B, 1-3) and the total synaptic current (I_syn-e; C, 1-3) recorded from each of the three cells in the simulation of Fig. 12 under a simulated voltage clamp at a holding potential of 60 mV. Because the bursting cells are coupled globally, the time courses of synaptic currents in B, 1-3, are similar. A comparison of the membrane potential trajectory with _syn-e for cell 1 reveals that I_syn-e has shown no appreciable activation when the cell has started firing. Thus in the case of the pre-I cells, intrinsic properties appear to determine the onset of spiking, and the subsequent increase in spike frequency when a majority of the population fires is due to a combination of both synaptic and intrinsic (voltage-dependent activation of I_NaP-h) mechanisms. In contrast in cells 2 and 3, the onset of the burst occurs after a noticeable increase in I_syn-e, suggesting a greater role for synaptic mechanisms in initiating the burst, especially in cell 3.

The follower cell populations were heterogeneous in E_L, and each cell also received input from different randomly selected neurons in the bursting population. Figure 13, D-F, illustrates the spike-frequency histograms (D, 1-3), typical membrane potential trajectory (E, 1-3), and a typical (gray) and average (black) I_syn-e elicited by a simulated voltage clamp at 60 mV for each