Models of Respiratory Rhythm Generation in the Pre-Bötzinger Complex. II. Populations of Coupled Pacemaker Neurons

Robert J. Butera Jr., John Rinzel, Jeffrey C. Smith


We have proposed models for the ionic basis of oscillatory bursting of respiratory pacemaker neurons in the pre-Bötzinger complex. In this paper, we investigate the frequency control and synchronization of these model neurons when coupled by excitatory amino-acid-mediated synapses and controlled by convergent synaptic inputs modeled as tonic excitation. Simulations of pairs of identical cells reveal that increasing tonic excitation increases the frequency of synchronous bursting, while increasing the strength of excitatory coupling between the neurons decreases the frequency of synchronous bursting. Low levels of coupling extend the range of values of tonic excitation where synchronous bursting is found. Simulations of a heterogeneous population of 50–500 bursting neurons reveal coupling effects similar to those found experimentally in vitro: coupling increases the mean burst duration and decreases the mean burst frequency. Burst synchronization occurred over a wide range of intrinsic frequencies (0.1–1 Hz) and even in populations where as few as 10% of the cells were intrinsically bursting. Weak coupling, extreme parameter heterogeneity, and low levels of depolarizing input could contribute to the desynchronization of the population and give rise to quasiperiodic states. The introduction of sparse coupling did not affect the burst synchrony, although it did make the interburst intervals more irregular from cycle to cycle. At a population level, both parameter heterogeneity and excitatory coupling synergistically combine to increase the dynamic input range: robust synchronous bursting persisted across a much greater range of parameter space (in terms of mean depolarizing input) than that of a single model cell. This extended dynamic range for the bursting cell population indicates that cellular heterogeneity is functionally advantageous. Our modeled system accounts for the range of intrinsic frequencies and spiking patterns of inspiratory (I) bursting cells found in the pre-Bötzinger complex in neonatal rat brain stem slices in vitro. There is a temporal dispersion in the spiking onset times of neurons in the population, predicted to be due to heterogeneity in intrinsic neuronal properties, with neurons starting to spike before (pre-I), with (I), or after (late-I) the onset of the population burst. Experimental tests for a number of the model’s predictions are proposed.


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).


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 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−6 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:

Run the simulation for 60 s (of simulation time) to allow transients to decay.
Collect interspike intervals (ISIs) for 60 s. Let pequal the maximum of these ISIs. Let w define a window size equal to 0.9p.
If there were no spikes, the cell is classified as silent.
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 sizew 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.
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.
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 asbeating. Otherwise, the cell is classified asbursting. 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 ofs 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.

Fig. 1.

Effect of level of depolarization (viag tonic-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. Forg tonic-e = 0.2 and syn-e = 15 nS the coupled system is bistable, with coexisting silent and beating solutions (top left). Parameter values of g tonic-eand 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 beCdVdt=INaP­hINaIKILItonic­eIsyn­e+Iapp Equation 1The 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(VE 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.

View this table:
Table 1.

Parameters for randomized normal distribution of intrinsic parameter values used in illustrated simulations

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 asIsyn,j=i=1Nḡsyn­e,i,jsi(VjEsyn­e,j) Equation 2where the gating variable s is driven by presynaptic depolarization according tos˙i=[(1si)s(Vi)krsi]/τs Equation 3This 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 functions (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 functions (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).

View this table:
Table 2.

Parameters for randomized normal distribution of synaptic conductances used in illustrated simulations


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 andg 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-ewere 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-eas g tonic-e is varied are plotted in Fig.2.

Fig. 2.

Effects of g tonic-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. 2where bursting activity occurs, several trends are evident.

At a given level of g tonic-e, syn-e increases both the burst duration (Fig. 2 A) and the burst interval (Fig. 2 B), 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 ofI 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 bySomers 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 withg 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 ofg 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 betweeng tonic-e and syn-ethat 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 withg 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 ofg 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.2 A), 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 Ldetermines 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 inButera 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 solelyE 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. Theleft column illustrates the dynamics of the uncoupled network ( syn-e = 0), and theright 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. 3 A1). There is little correlated network activity (population spike histogram, Fig. 3 A2). This dispersion in intrinsic activities is consistent across 50 simulations with different randomly generated parameter sets and initial conditions (bar plots, Fig.3 A3). After coupling, most of the cells burst in relative synchrony with a coordinated periodic burst of action potentials (Fig.3 B, 1 and 2). This finding is also consistent across a large number of simulations (B3).

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 Ncycles, 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. 4 C).

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-Ca2+ 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-Ca2+solution (uncoupled), allowing an assessment of the role of coupling on the dynamics of bursting at a cellular level. Figure5 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).

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 Table3. 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 (Tables4 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 this table:
Table 3.

Parameter set distributions used for various coupling simulations

View this table:
Table 4.

Distribution of activity modes in uncoupled and coupled heterogeneous populations

View this table:
Table 5.

Effects of synaptic coupling on burst frequency and burst duration

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 ofE 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 7 A, 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.

Fig. 7.

Voltage-dependent frequency control of network activity.Left: histograms of network spiking (10-ms bin size) as the mean level of g tonic-e is increased.B: network burst frequency as a function ofg tonic-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. 7 A, 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.7 A7), 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 ofE L and NaP-h and identical g tonic-e values is illustrated in Fig.7 C. This reveals that the dynamic range ofg 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 theg 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 ofg tonic-e larger than the single cell but less than that for the heterogeneous network (Fig. 8 A). 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. 8 B). From these results, we conclude that both intrinsic heterogeneity and excitatory synaptic connectivity contribute toward the increased dynamic range ofg 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.

Fig. 8.

Effects of parameter heterogeneity and coupling strength on frequency control of network bursting. Curves plot burst frequency as a function of g tonic-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 (E L, NaP-h, syn-e) with networks with intrinsic heterogeneity only (E L, 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. 8 C). 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 bothg 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 ofE 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 9also 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 (seediscussion).

Fig. 9.

Intrinsic modes of activity for disconnected component cells in simulation of Fig. 7. A: each horizontal line represents the range of g tonic-e where individual cells exhibit intrinsic bursting. B: number of cells demonstrating intrinsic silence, bursting, or beating activity modes as a function of g tonic-e. C: ◊, mean, minimum, and maximum burst frequency of cells with intrinsic bursting activity as a function ofg tonic-e. For comparison, ⧫, synchronous burst frequency as a function ofg tonic-e for the coupled network.

Emergent network oscillations from coupled cells with INaP-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 orI 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 Ldetermining 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. 10 A, the uncoupled population shows only five neurons exhibiting spiking activity, the rest are silent. Coupling the population (Fig. 10 B) 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. 10 C) across the entire population. Alternatively, an increase in the mean value of NaP-h (Fig. 10 E) 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. 10 D). 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).

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 E L orI app. 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. 8 C(×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:

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.
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 withNaP-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 ELyields 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.
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 withNaP-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 ELand syn-e are specified in Tables 1 and2. 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.11 A, 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.11 B, 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. 11 A 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 anidentical assignment of intrinsic and synaptic parameters, only using a mean synaptic conductance (0.19 nS) viag 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).

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. 11 C, 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. 11 C. 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-Ca2+ 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.

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.11 C 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.

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. 12 B 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 inB, 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 ofI 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 inE L, and each cell also received input from different randomly selected neurons in the bursting population. Figure13, 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 cell (F, 1–3).

These cells have no intrinsic bursting properties and fire action potentials when synaptic input exceeds a given amount. The follower cells simply “read-out” the dynamics of the population activity of bursting cells providing their synaptic input, as reflected in the time course of the synaptic drive current envelope (Fig.13 F). The spike-frequency histograms of these follower cells have a similar spectrum of shapes, ranging from pre-I (D1–F1) and I (D2–F2) to late-I (D3–F3), as those of the bursting cells in Fig. 12. These spiking patterns have been found (Johnson et al. 1994) for neurons without intrinsic bursting properties in the pre-BötC of neonatal rat slice preparations. The rapidly peaking, slowly decaying time course of the synaptic current envelope (Fig.13 F) is characteristic of follower-type neurons in the slices (Funk et al. 1993; Smith et al. 1991) and the more en bloc in vitro preparations (Smith et al. 1990, 1992). Our results suggest that a range of both intrinsic (for bursters) and synaptic (for bursters and followers) heterogeneous properties give rise to the spectrum of spike-frequency histograms recorded from pre-BötC inspiratory neurons.

The amplitude of I syn-e (Fig. 13 C) in bursting cells is smaller than that required to generate subthreshold depolarization and the firing of action potentials in the follower cells (Fig. 13 F). The only difference between the follower cells and the bursting cells is the absence ofI NaP-h. The presence ofI NaP-h not only allows intrinsic bursting to occur but means that less additional depolarizing input (such asI syn-e) is required to elicit a burst. This postulated role of I NaP-h has been referred to as synaptic amplification (e.g., Ramirez and Richter 1996), which has been hypothesized (Ramirez and Richter 1996; Rekling and Feldman 1998; Smith et al. 1995) to promote synchronized bursting activity in the pacemaker population. Our simulations demonstrate this function ofI NaP-h, which is particularly apparent in cases where most of the burst-capable cells in the population are silent but fully synchronous bursting occurs readily with synaptic input (e.g., Figs. 9 and 10).


In this study, we have investigated the dynamics of populations of model pacemaker cells incorporating postulated subthreshold mechanisms for intrinsic burst generation and excitatory synaptic mechanisms for the synchronization and control of respiratory rhythm. We have hypothesized (Smith 1997; Smith et al. 1995) that this hybrid of pacemaker and network mechanisms is the basis of rhythm generation in the pre-BötC. It has not been proven that a population of the type of cells we have modeled is the actual rhythm-generating kernel in the pre-BötC although the mechanisms we have incorporated in the model are consistent with in vitro experimental data (reviewed in Rekling and Feldman 1998; Smith et al. 1995). Experimental data on intrinsic currents and network architecture is lacking. Our approach therefore has been to explore the dynamics of populations of minimal models of cells with plausible burst-generating mechanism(s) and to investigate the consequences of heterogeneous cellular properties and connectivity schemes incorporating heterogeneous distributions of connection strengths as well as varying degrees of sparsity of connections for small and large populations of cells. Our objective was to determine general principles of network operation that may account for available experimental observations on patterns of bursting pacemaker and follower cell activity in the network. These general principles include:

excitatory synaptic tonic drive increases the synchronous bursting frequency of the population, while phasic excitatory synaptic coupling decreases the frequency.
in a population of burst-capable cells (some degree ofNaP-h is present in all cells), only a minority (as few as 10%) of the population needs to demonstrate intrinsically bursting or beating activity modes to trigger network-wide synchronized oscillatory bursting; in cases of lowNaP-h, synchronized oscillations can arise as an emergent property of the network.
sparsity of excitatory coupling (as low as 3%) has little effect on the synchronization of individual network-wide bursting, however, sparse networks exhibit greater variability in burst-to-burst cycle length.
except at very large values of syn-e, the population of cells exhibits synchronized bursting over a wider dynamic range of input parameters (such as gtonic-e) than that of a single cell, and both excitatory coupling and intrinsic heterogeneity contribute to this effect.
increasing the strength of excitatory coupling (for values ofsyn-e > 1 nS) served to decrease the range of gtonic-e where bursting occurred, as well as the range of synchronous burst frequencies.
asynchronous tonic input from a finite population of beating cells, as opposed to a mean level of tonic input, increases the number of cells that burst when uncoupled.
both burster and follower cells may exhibit the range of spike-frequency histograms (pre-I, I, late-I) recorded from inspiratory neurons in pre-BötC slice preparations.
bursting neurons with higher intrinsic burst frequencies tend to be more depolarized and fire earlier in the cycle than bursting neurons with lower intrinsic frequencies.

Comparison with experimental data and model limitations

Johnson et al. (1994) studied the effect of synaptic coupling on pre-BötC bursting neuron activity by examining burster cell dynamics with synaptic transmission intact and blocked by low-Ca2+ solution. We modeled these experiments by analyzing the dynamics of a heterogeneous population of model neurons before and after elimination of EAA-mediated synaptic coupling. Our model predictions are consistent with experimental observations. However, there are several possible flaws with the technique applied inJohnson et al. (1994) and our approach in modeling these experiments. First, the use of a low-Ca2+ solution to eliminate synaptic connectivity makes it possible that Johnson et al. also observed the effect of altering the dynamics of Ca2+currents in the bursting neurons they were measuring. Similar experiments have been performed using 6-cyano-7-nitroquinoxaline-2,3-dione (CNQX) (Smith et al. 1997) to block EAA-mediated transmission, and similar quantitative effects (Figs. 5 and 6) on burst dynamics have been observed. Second, blocking synaptic transmission eliminates not only synaptic coupling between bursting neurons but also synaptic input from tonic cells onto the population of bursters. We have not modeled this effect because the amount of excitatory tonic input received by a bursting pre-BötC cell has not been quantified. The effect of excitatory tonic input on single model cells is to depolarize them and increase their frequency (see Fig. 1) (see also Fig. 6 of Butera et al. 1999), whereas the effects of EAA coupling between bursting neurons is to decrease their frequency. Figure 9 illustrates that simultaneous removal of synaptic coupling within the population as well as a modest reduction in depolarizing input (modeled byg tonic-e) still will result in an increase in the mean bursting frequency of the uncoupled population. Third, we have not explicitly modeled inhibitory synaptic connections from a population of tonic inhibitory cells that also may provide modulatory input in vitro (see discussion in Smith et al. 1995). The experimental data and our model predictions show that synaptic connectivity leads to a reduction in mean bursting frequency and increase in burst duration; thus we conclude that the effect of EAA-mediated coupling among bursting neurons is more significant (under the experimental conditions of Johnson et al. 1994) than the effect of EAA-mediated or inhibitory tonic drive to the bursting neurons. Fourth, bidirectional electronic coupling has been reported in respiratory motoneurons (Rekling and Feldman 1997). Although it has not been identified in pre-BötC interneurons, it remains a potential coupling mechanism. However, Koshiya and Smith (1998, 1999) have shown using Ca2+imaging that pre-BötC neurons with pacemaker properties continue to burst, but desynchronize, in the presence of CNQX. Thus even if electrotonic coupling is present, it may not be a major mechanism for synchronizing the activity within the slice. Electrotonic coupling is also unlikely to account for bilateral synchronization of the rhythm.

Our model neuron populations have temporal patterns of neuronal spiking as shown by spike-frequency histograms, as well as patterns of synaptic drive currents in individual cells, similar to those found experimentally for burster and nonburster follower neurons in vitro. However, our neurons are simplified one-compartment models designed to examine how postulated subthreshold mechanisms for intrinsic burst generation can account for the dynamics of both individual neurons and populations of bursting neurons in the pre-BötC. We have deliberately not attempted to model the role of other specific ionic currents (such as I A, I T, or I K,Ca) in shaping the firing patterns of action potentials in pre-BötC neurons, although they certainly can modulate the firing patterns of isolated respiratory neurons as well as those embedded in a network (e.g., Rybak et al. 1997a,b). Interestingly, in the case of follower-type cells in vitro, we previously have shown that the shape of spike frequency histograms and time course of synaptic drive potentials largely reflect the time course of the excitatory synaptic drive current envelope (e.g., Funk et al. 1993; Liu et al. 1990). The time course of excitatory drive current (rapidly peaking, slowly decaying) exhibited by follower cells in our model is similar to that seen experimentally and is reflective of the overall time course of the presynaptic cell population activity which is read out by the follower cells. We have also not included in our model dendritic compartments even though respiratory neurons may have a significant number of proximal dendrites (e.g., Funk et al. 1993), which could affect the ionic current distribution and firing patterns in response to synaptic drive. More complex models that include the effects of ionic currents that modulate firing frequency and dendritic compartments will provide addition insight into the shape and amplitude of the spike-frequency histograms. More quantitative information is required from experimental measurements on the types and kinetics of ionic currents as well as morphology intrinsic to respiratory pre-BötC neurons.

Although our single cell models possess membrane capacitance values and resting membrane impedances consistent with whole cell measurements from respiratory neurons in vitro (Smith et al. 1991,1992), the magnitude of the synaptic or applied currents required to change membrane potential of the bursting cells may be somewhat smaller in the model than typically observed experimentally (N. Koshiya and J. Smith, unpublished observations). We attribute this mainly to I NaP-h in the model, whose voltage-dependencies and maximal conductances have not been determined, and these factors greatly affect the amount of current input required to elicit a burst.

Model predictions and interpretations


This study has illustrated that both excitatory synaptic coupling and intrinsic parameter heterogeneity make rhythmic bursting more robust. Figure 2 illustrates that with identical cells, EAA-mediated coupling increases the range of g tonic-e where synchronous bursting is supported as opposed to the case of a single cell ( syn-e = 0). Figures 7 and 8 show similar effects when comparing the range ofg tonic-e where bursting is supported in populations as compared with single cells. The functional consequences of heterogeneity in bursting neuron networks have not been analyzed previously. Smolen et al. (1993) studied models of clusters of heterogeneous populations of pancreatic β-cells. While they considered a different form of connectivity (3-dimensional, local, gap-junctional coupling), a similar finding was made: the range of a key control parameter (a glucose sensing conductance) where bursting occurred was larger for the coupled population than for individual cells. Thus it is plausible that the inherent cell-to-cell variability in intrinsic properties is a functionally advantageous factor in maintaining robust rhythmic output from a network of cells with intrinsic bursting properties. Information on the degree of heterogeneity in the pre-BötC pacemaker cell network remains to be obtained. We predict that in vitro, the range of values of a depolarization parameter (e.g., [K+]o) over which bursting occurs is broader for the respiratory network in slices containing the pre-BötC than for individual isolated inspiratory bursting neurons within the pre-BötC.


Our model simulations have illustrated the principle of voltage-dependent frequency control of both isolated inspiratory pre-BötC pacemaker neurons (Butera et al. 1999) and an excitatory-coupled network of modeled pre-BötC neurons (Fig. 7). In both cases, the frequency of the bursting oscillation is regulated by the amount of depolarizing input, such asg tonic-e or E L. Similar effects have been shown at the level of a single cell viaI app (Smith et al. 1991) or perturbations in [K+]o (R. Butera, C. Del Negro, and J. Smith, unpublished data). We predict that titrating the level of excitation and frequency of bursting in the pre-BötC slice (e.g., by controlling [K+]o) will show the following effects: an increase in frequency with elevated [K+]o, a decrease in overall amplitude of burster cell and thus follower cell population activity, and a change in the shape of the integrated population activity from a rapidly peaking population burst at hyperpolarized potentials to a more gradual incrementing-decrementing profile. Assuming that the integrated XII motor output is an accurate measure of overall network spiking activity in the slice, changes in XII motor output similar to the spike-activity plots of Fig. 7 A, 1–7, should be obtained as pre-BötC neuron tonic excitation is increased.


The simulations show that with a heterogeneous pacemaker cell population, many of the neurons may not be spontaneously bursting when isolated but are burst-capable, and these burst-capable cells may be transformed into phasically bursting neurons in the presence of phasic synaptic input. The firing of these cells further contributes to the robustness of population bursting by generating additional excitatory synaptic input to the population.

We predict that, in the extreme case, it is not necessary for any neurons to exhibit intrinsic oscillatory bursting for synchronous bursting to emerge in the network (Fig. 10). This possibility also has been suggested by Rekling and Feldman (1998), although not rigorously analyzed. What our simulations indicate is necessary is that a significant fraction of the neurons in the network possess some degree of I NaP-h (or some other subthreshold cationic current capable of synaptic amplification, see following text) and that a small fraction (≥10%) must exhibit intrinsic neuronal activity whether that activity be spiking or bursting. It is currently unknown whether these conditions occur physiologically. We assume that the network normally operates in vitro in a mode with some cells in the pacemaker population in their oscillatory bursting state. The experimental data from slice preparations is consistent with this assumption, at least for the conditions of the slice, because some pre-BötC inspiratory neurons continue to exhibit rhythmic bursting after block of synaptic transmission (Johnson et al. 1994; Smith et al. 1997). We previously have assumed that many of the bursting neurons would be in their silent state (Smith 1997; Smith et al. 1995), but we currently do not know the fraction of the population that is sufficiently hyperpolarized or has sufficiently low NaP-h to be in a nonoscillatory state and must be recruited into synchronized bursting by synaptic input. We have proposed several in vitro tests (Butera et al. 1999) that may identify burst-capable cells, e.g., a poststimulus burst may be triggered by a hyperpolarized voltage or current clamp or a tonic cell illustrates spike-frequency adaptation (initially high rebound firing frequency) after a prolonged hyperpolarization.

Our modeling results also predict that pacemaker type pre-BötC neurons would require less synaptic input current to induce bursting, compared with follower-type pre-BötC neurons that do not possessI NaP (Fig. 13). This is consistent with preliminary experimental observations (Smith et al. 1991; N. Koshiya and J. Smith, unpublished observations) and may be explained by the notion of synaptic amplification (Crill 1996; Ramirez and Richter 1996; Schwindt and Crill 1995), where the presence of a subthreshold inward current like I NaP means that less depolarizing synaptic input is required to initiate a burst.


Synchronized bursting is predicted to occur with substantial heterogeneity, but desynchronization can occur with very weak coupling, extreme parameter heterogeneity, and/or low levels of depolarizing input. Recent experiments on pre-BötC slice preparations in mice (Ramirez et al. 1996, 1998) have reported N:1 (such as 3:1) coupling patterns between pre-BötC neurons and hypoglossal motor output in slices from juvenile mice. We predict that a breakdown of 1:1 coupling would occur in the pre-BötC slice preparations when synaptic coupling strength is progressively reduced (e.g., by CNQX). Such a breakdown in synchronized bursting may be tested for by using activity-dependent imaging techniques to visualize bursting activity of multiple pacemaker cells simultaneously (Koshiya and Smith 1998, 1999).


We put forward the notion that the different spike-frequency histograms recorded from oscillatory bursting neurons in pre-BötC slice preparations may not be due to differences in the existence of specific types of ion channels or synaptic connectivity but rather may be accounted for by a single model with variable levels of expression of ionic currents responsible for generating bursting dynamics. Pre-I, I, and late-I cells simply may possess different levels of expression of ion channels and/or synaptic input regulating intrinsic levels of depolarization and the ability for burst initiation. The possible role of pre-BötC neurons with pre-I spiking patterns in rhythm generation has been discussed (e.g., Schwarzacher et al. 1995; Smith et al. 1997). The early low frequency spiking of the pre-I cells provides excitation and contributes toward recruiting the bursts of the remaining cells in the network; they do not dictate the overall frequency of the network, however, because many other cells in the network are also intrinsic pacemakers or burst-capable and collectively set the oscillation frequency. Thus the pre-I cells cannot be considered as simple triggers of the I-phase. The frequency of the coupled network is (in all the simulations we have studied) lower than the frequency of the pre-I cells when uncoupled. We predict that pre-I neurons should correspond to neurons with higher intrinsic bursting frequency or baseline levels of depolarization (or both), while neurons which fire later in the cycle should be more hyperpolarized and/or possess lower intrinsic frequencies.


We have analyzed the dynamics of excitatory networks of voltage-dependent pacemaker neurons as a model for the rhythm and inspiratory burst generating kernel in the pre-BötC of in vitro preparations. The general principles of operation of this pacemaker network have been investigated. We conclude that a heterogeneous network of the type of voltage-dependent burster cells modeled can account for the synchronization, frequency control, patterns of inspiratory cell activity, synaptic interactions, and network activity found in vitro. The voltage-dependent properties of the bursting cells combined with tonic drive and phasic excitatory interactions can provide a functionally robust mechanism for frequency control over a wide dynamic range. The model requires further experimental verification, and we have presented predictions to facilitate additional tests. At present the model is intended to account for patterns of inspiratory activity generated in vitro. An important extension of the model will be to embed the pacemaker kernel in a pattern formation network (Smith 1997) that provides additional synaptic control, including phasic inhibition. Analyzing dynamic interactions between rhythm and pattern generating elements that produce the complete pattern of inspiratory and expiratory phase activity in vitro and in vivo remains an important problem.


We thank N. Koshiya, S. Johnson, C. Wilson, and G. de Vries for helpful discussions and R. Burke and A. Sherman for critical readings of the manuscript. J. Rinzel thanks the Laboratory of Neural Control, NINDS, for hosting him in the summer of 1998.

This work was supported by the intramural research programs of National Institute of Neurological Disorders and Stroke and National Institute of Diabetes and Digestive and Kidney Diseases.

Present address of R. J. Butera, Jr.: School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA 30332-0250.


  • Address for reprint requests: J. C. Smith, Laboratory of Neural Control, NINDS, NIH, Bldg. 49, Room 3A50, 49 Convent Dr., Bethesda, MD 20892-4455.

  • 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.


View Abstract