Intrinsic Bursters Increase the Robustness of Rhythm Generation in an Excitatory Network

L. K. Purvis, J. C. Smith, H. Koizumi, R. J. Butera


The pre-Botzinger complex (pBC) is a vital subcircuit of the respiratory central pattern generator. Although the existence of neurons with pacemaker-like bursting properties in this network is not questioned, their role in network rhythmogenesis is unresolved. Modeling is ideally suited to address this debate because of the ease with which biophysical parameters of individual cells and network architecture can be manipulated. We modeled the parameter variability of experimental data from pBC bursting pacemaker and nonpacemaker neurons using a modified version of our previously developed pBC neuron and network models. To investigate the role of pacemakers in networkwide rhythmogenesis, we simulated networks of these neurons and varied the fraction of the population made up of pacemakers. For each number of pacemaker neurons, we varied the amount of tonic drive to the network and measured the frequency of synchronous networkwide bursting produced. Both excitatory networks with all-to-all coupling and sparsely connected networks were explored for several levels of synaptic coupling strength. Networks containing only nonpacemakers were able to produce networkwide bursting, but with a low probability of bursting and low input and output ranges. Our results indicate that inclusion of pacemakers in an excitatory network increases robustness of the network by more than tripling the input and output ranges compared with networks containing no pacemakers. The largest increase in dynamic range occurs when the number of pacemakers in the network is greater than 20% of the population. Experimental tests of our model predictions are proposed.


Networks of neurons that generate rhythmic networkwide bursts of action potentials are found in many areas of the nervous system. One such network is the pre-Botzinger complex (pBC), a population of neurons that is critical for generating the respiratory rhythm (Gray et al. 2001; Smith et al. 1991). Despite extensive research to date, the details of the mechanisms of rhythmogenesis in the pBC are still not completely resolved. Many experiments have established that the pBC both in vitro and in more intact states contains intrinsic bursters or “pacemaker” neurons—cells that are capable of rhythmic bursting in the absence of synaptic input (Del Negro et al. 2002a, 2005; Johnson et al. 1994; Koshiya and Smith 1999; Pagliardini et al. 2005; Paton et al. 2006; Pena et al. 2004; Smith et al. 1991; Thoby-Brisson and Ramirez 2001; Tryba et al. 2006). Two types of pacemaker neurons are found in the pBC: one is dependent on a persistent sodium current (NaP dependent, also cadmium insensitive; Del Negro et al. 2002a,b, 2005; Paton et al. 2006; Pena et al. 2004; Thoby-Brisson and Ramirez 2001; Tryba et al. 2006) and the other uses a calcium-dependent mechanism (CaN dependent, also cadmium sensitive; Del Negro et al. 2005; Pena et al. 2004; Thoby-Brisson and Ramirez 2001; Tryba et al. 2006). Furthermore, these and other pBC neurons are coupled by excitatory synaptic connections (Funk et al. 1993; Koshiya and Smith 1999) and the pBC excitatory network itself exhibits autorhythmic properties (Johnson et al. 2001). Recent literature raises questions about the significance and the abundance of pacemaker cells within the pBC network (Del Negro et al. 2002a,b, 2005; Pagliardini et al. 2005; Paton et al. 2006; Pena et al. 2004; Tryba et al. 2006). In the intact respiratory network, the pBC is part of a much larger circuit (Feldman and Del Negro 2006; Smith et al. 2000) and pBC cells are responsive to modulation that controls inspiratory frequency (Johnson et al. 1996; Schwarzacher et al. 2002; Solomon et al. 2000, 2002). It was also established that the number of pacemakers in the pBC can be dynamically varied by neuromodulators acting on the network (Pena and Ramirez 2002, 2004; Ramirez et al. 2004). Here we use modeling studies to investigate how the fraction of the network population that is composed of intrinsic bursters affects measures of rhythmogenic capacity or “robustness” of network activity. Robustness is quantified in terms of both the input range and output range of the network. The input range is the size of parameter space where networkwide synchronous rhythmic bursting occurs and the output range is the range of bursting frequencies that the network produces across this input range. Thus a network is considered robust if it bursts at many levels of tonic excitation and over a wide range of frequencies. In general, the question of the contributions of pacemaker neurons to the dynamic performance of excitatory networks was not previously addressed quantitatively.

We ran network simulations to investigate this problem using a modified version of our previously developed pBC neuron and network models (Butera et al. 1999a,b). We simulated a heterogeneous excitatory network consisting of 50 pBC neurons—the smallest number of neurons that captures the dynamics of larger networks (Butera et al. 1999b)—and we varied the fraction of the population made up of pacemakers from 0 to 100%. While varying the fraction of pacemakers, we also varied the amount of excitatory drive supplied to the network at several different levels of excitatory synaptic coupling strength. Simulations with no pacemakers in the network were able to produce synchronous networkwide rhythmic activity; however, the input and output ranges were substantially reduced. Our simulations suggest that the inclusion of pacemaker neurons in a bursting network allows for greater controllability of the frequency range produced by that network. We conclude by proposing experimental tests of our model predictions.


Experimental data

Experimental data obtained from neonatal rat pBC inspiratory neurons in vitro shows from voltage-clamp recording analysis that a primary difference between pacemakers (PMs) and nonpacemakers (NPMs) in the pBC is the ratio of the persistent sodium conductance to the leak conductance (gNaP/gLeak): to date all pBC inspiratory neurons analyzed in vitro exhibit NaP and K+-dominated leak currents, although the gNaP/gLeak ratio is typically larger in PMs than in NPMs (Del Negro et al. 2002a; Koizumi and Smith 2004). When viewing the recent data of Koizumi and Smith (2004) as a scatterplot in gNaP and gLeak space (Fig. 1 A), there is a clear distinction between the PMs and NPMs (dashed line). Because the hand-drawn dashed line in Fig. 1A has a nonzero y-intercept, the gNaP/gLeak ratio is not a perfect indicator of the single-cell dynamics (i.e., whether the cell is a PM or NPM). To accurately compare the experimental data to the model (see following text) in terms of conductance densities, the data for individual neurons (n = 50) were normalized by the measured capacitance and then multiplied by the nominal whole cell capacitance used in the model (21 pF). PMs used for this analysis were differentiated from NPMs after pharmacological block of glutamatergic excitatory synaptic currents as previously described (Del Negro et al. 2002a; Koshiya and Smith 1999) and the dependency of intrinsic rhythmic bursting on gNaP was confirmed pharmacologically for identified PMs. All gNaP and gLeak data used in this study are from the set of PMs and NPMs obtained by Koizumi and Smith (2004) and are reproduced here for data parameter space mapping, model parameter adjustments, and data-model comparisons.

FIG. 1.

Parameter space and pacemaker/nonpacemaker (PM/NPM) boundaries for experimental data and the model. Pacemakers are represented by, nonpacemakers by *. A: scatterplot of experimental data (Koizumi and Smith 2004) obtained from voltage-clamp analysis of PM/NPM neurons identified in the pre-Botzinger complex (pBC) of neonatal rat in vitro slice preparations, using neuron visualization and recording techniques similar to those described previously (Del Negro et al. 2002a; Koshiya and Smith 1999). Note the clear boundary between pacemakers and nonpacemakers (hand-drawn dashed line). n = 30 PMs; n = 20 NPMs. B: model's PM/NPM boundaries. Parameter set was classified as a PM if the model neuron bursts at any level of stimulus current (dark region).

Neuron model

This study uses our previous pBC neuron model (see appendix) and heterogeneous excitatory network models of these cells (Butera et al. 1999a,b). The neuron model (Model 1 of Butera et al. 1999a) is a Hodgkin–Huxley (Hodgkin and Huxley 1952) style, electrophysiological model of a bursting pBC cell. The model consists of three voltage-gated ionic currents and a K+-dominated ohmic leak current. The voltage-gated currents are 1) a fast Na+ current, 2) a delayed-rectifier K+ current, and 3) a slowly inactivating persistent Na+ current (NaP). NaP is responsible for the intrinsic voltage-dependent, rhythmic bursting behavior displayed by these neurons (i.e., this is a model of the NaP-dependent or cadmium-insensitive bursters; Butera et al. 1999a). This is the biophysically minimal set of currents required to describe the main features (below) of pBC neuron properties. Other currents include a tonic excitatory drive current and an excitatory synaptic current. The excitatory drive, modeled as gtonic, models input from a tonic spiking population, which is proposed to function for pBC network burst frequency control (e.g., see Butera et al. 1999b; Smith et al. 2000). The excitatory synaptic current models the excitatory amino acid–mediated coupling between individual bursting neurons in all simulations using the conductance gsyn (Butera et al. 1999b).

Some parameters of the original Model 1 neurons were adjusted to more closely match recent experimental data. Specifically, more accurate values for NaP half-activation (V1/2max = −45 mV) and slope factor (k = 5) are used (Koizumi and Smith 2004). Because the NaP half-activation value was hyperpolarized by 5 mV, we also hyperpolarized the NaP half-inactivation by 5 mV. These adjustments are also consistent with other measurements of NaP properties of neurons isolated from the pBC region in vitro (Ptak et al. 2005; Rybak et al. 2004). NaP inactivation kinetics was not previously quantified experimentally for pBC neurons, so the original Model 1 kinetics was used. Simulations with these adjusted models provide behavior very similar to that of the original Model 1 neurons including voltage-dependent rhythmic bursting with a similar range of oscillation frequencies, controllable by applied current or tonic excitatory synaptic input (see Fig. 6 in Butera et al. 1999a). Figure 2 illustrates the pacemaker and nonpacemaker behaviors exhibited by the model for two different values of gNaP. NPMs make the transition from silence directly to beating as excitatory drive is increased (Fig. 2A, gNaP = 1.5 nS), whereas PMs make the transition from silence to bursting to beating (Fig. 2B, gNaP = 2.5 nS). This model also captured prominent properties of the recorded data neurons, justifying our parameter sets, including 1) silent, rhythmic bursting, and tonic spiking regimes determined by baseline membrane potential as controlled by an applied current; 2) bursting frequencies tunable over an order of magnitude range of frequencies by applied current (see also Koshiya and Smith 1999); 3) subthreshold current–voltage relations obtained during slow voltage-clamp ramps (30 mV/s used to obtain the data) that are separable into gNaP and gLeak as the two main conductance components (see also Butera et al. 2005); and 4) rhythmic bursting controllable as seen experimentally by gNaP and gLeak (i.e., by gNaP/gLeak ratios). Detailed comparisons of the data-model gNaPgLeak parameter spaces (Fig. 3) also indicate that our simplified model is consistent with the differentiation between data PMs and NPMs in the majority of cases.

FIG. 2.

Model behavior for 2 different values of gNaP. A: gNaP = 1.5 nS. NPMs progress from silence to beating (tonic firing) as gtonic is increased. B: gNaP = 2.5 nS. PMs progress from silence to bursting to beating as gtonic is increased. Only gNaP is changed (gLeak = 2.2 nS in A and B) to transform the NPM in A into the PM in B. *Transient firing activity in A is caused by step increase of gtonic and is not a burst.

FIG. 3.

Comparing data space to model space. PMs represented by, NPMs by *. Experimental data from Fig. 1A are corrected by increasing gNaP values by 25%. A: corrected experimental data plotted with the model's PM/NPM boundary. B and C: plots of data and random parameter distributions with the outline of the conservative PM region (denoted by solid lines) and NPM region (denoted by dashed lines) of model parameter space. In B, experimental data are plotted with the conservative model boundaries. Several data points fall outside of the conservative boundary allowed for model parameters. C: illustration of a sample randomly generated model parameter distribution with 30 PMs and 20 NPMs.

The network simulations required a mixed population of models of PMs and NPMs. Although previous modeling results classified the activity modes of the model (beating, bursting, silence; see Fig. 7 of Butera et al., 1999a) at rest (i.e., with no externally applied current), the PM and NPM classification of Fig. 1A was experimentally determined by testing for bursting at any level of applied current under current-clamp recording. Therefore a similar process was applied to the model using the interactive differential equation simulation package XPP (Ermentrout 2002) for model simulations along with Perl scripts and MATLAB (The MathWorks, Natick, MA) for data analysis. Both gNaP and gLeak were swept from 0 to 6 nS, which covers the full range of experimentally measured values. For each set of parameters represented in Fig. 1B, the stimulus current was swept from −30 to +30 pA. The stimulus current was started at −30 pA to ensure that the cell is at rest. If the values of gNaP and gLeak are such that the model is a nonpacemaker, then the modes of activity as the stimulus current is increased will progress from silence directly to beating (Fig. 2A). If rhythmic bursting occurs during part of this range (Fig. 2B), the neuron is classified as a pacemaker with this parameter set. The results of these simulations are given in Fig. 1B. Figure 1, A and B shares a clear boundary between the PM and the NPM regions. When creating a population of PMs and NPMs for the network simulations, the parameters must be chosen based on the modes of activity defined in Fig. 1B.

Modeling intrinsic parameter variability

The heterogeneous network simulations use parameters chosen from a two-dimensional parameter space (only gNaP and gLeak are varied). For these simulations, we created two heterogeneous populations of PMs and NPMs whose variability reflects the variability that exists in the experimental data (Table 1). That is, we constructed population distributions reflecting not only the ranges of NaP and leak conductance densities but also the experimentally determined variability of these parameters (Table 1). Accomplishing this requires completing three tasks: 1) correcting experimental conductance data for measurement errors, 2) determining the model's acceptable operating range for PMs and NPMs, and 3) defining parameter distributions inside the model's acceptable operating range that provide mean and SD values consistent with the corrected experimental data.

View this table:

Conductance values for PMs and NPMs


It is known that experimental measurements of gNaP depend on the rate of voltage-clamp ramps (Del Negro et al. 2002a) used to estimate values of gNaP, which is consistent with our model kinetics for NaP inactivation. We simulated the ramp protocol used to determine the gNaP values of the experimental data we are attempting to model to produce an estimate of the underestimation resulting from the ramp rate. The simulation suggests an underestimation arising from the ramp rate of roughly 20% and we conservatively estimated an additional 5% error from other potential measurement errors (e.g., space-clamp). Therefore we corrected the gNaP data by multiplying by the expected underestimation of 25%. The value of gLeak was not found experimentally to be affected by ramp rates as predicted by the model and was therefore left unchanged. After correcting the data, 86% of the data points are on the appropriate side of the model's PM/NPM boundary (Fig. 3A).


To completely describe an allowable parameter space of the model for network simulations, both the PM/NPM boundary and the boundary that defines the operating range of the model must be defined. An accurate description of this space is necessary to reliably separate PM and NPM model neurons; parameter sets outside of the model's operating range will produce spurious results (e.g., plateau potentials) or nonphysiological models (e.g., negative conductances). The boundary between PMs and NPMs was chosen by fitting a diagonal line to the boundary previously computed (Figs. 1B and 3A). A model whose parameters lie on the PM/NPM boundary will produce doublet or triplet bursts of spikes. Therefore the y-intercept of the line was increased by 0.2 for PMs and decreased by 0.2 for NPMs (see white space between PM and NPM regions in Fig. 3C), to provide a conservative separation between PMs and NPMs and definitively determine that a model in that region is a PM or NPM. The minimum value for gNaP for all simulations was chosen to be 0.5 nS. The other diagonal boundary at the upper edge of the PM space was set to keep the model inside its operating range. Figure 3B displays these conservative model boundaries plotted with the experimental data of Fig. 1A. Seventy-four percent of the data points are inside the conservative boundaries allowed for the model parameters.


As stated earlier, the network simulations will use values from a two-dimensional normal distribution in gNaP and gLeak for PMs and NPMs. When randomly choosing parameters to generate a PM model, the parameters must be verified to ensure that the model actually is a PM (likewise for NPMs). The parameters must also be tested to ensure that they are not outside the operating range of the model. Thus some of the values that are randomly chosen from the two-dimensional distribution must be discarded because they fall outside these boundaries. Therefore the values for mean and SD of the experimental data listed in Table 1 cannot be used for the random-number generator because the act of discarding and repicking from the normal distribution will result in a different mean and SD for the population. C++ simulations were run that selected 10,000 random sets of parameters to determine the actual mean and SD values after bounding the allowable parameter space. The simulation was then repeated for various nominal means and SDs until the actual mean and SD after discarding/repicking was similar to the mean and SD of the corrected data (Table 1). These trials were repeated 100,000 times and the best parameter distributions were determined by taking the minimum of the summed percentage error between the computed mean/SD and the data mean/SD. Table 1 provides the optimal means and SDs for gNaP and gLeak whole distribution statistics that best match the experimental data.

Once initial mean/SD values for the normal distribution were chosen such that the calculated mean/SD after discarding/repicking closely matched the experimental data, we used these distributions to choose model parameters of PMs and NPMs for the network simulations. For example, Fig. 3C provides a sample parameter distribution to be used for model simulations that includes 30 PMs and 20 NPMs. Although this data-mapping scheme is not perfect (26% of the corrected experimental data falls outside the model's allowable range; Fig. 3B), it does provide a model approximation of the variability of the experimental data.

Network simulations

The network simulations consist of a population of 50 cells, of which there are 1K PMs and 50K NPMs. It is important to note that the total number of cells in the network remains constant at 50 while K is varied from 0 to 50 in increments of one (i.e., only the percentage of the population that are PMs or NPMs changes). For each K, the amount of excitatory drive, or gtonic, is varied in 0.1-nS increments from 0 to 1.5 nS. Varying the amount of drive allows us to explore this voltage-dependent frequency control, as exhibited by the isolated pBC in slice preparations (Del Negro et al. 2001; Smith et al. 1991). Also, because the exact amount of excitatory synaptic coupling is unknown, the level of excitatory coupling is varied over physiological ranges for each simulation. Physiological range was determined by examining the synaptic current for various values of gsyn (see, e.g., Fig. 13 of Butera et al. 1999b) and comparing it to values measured experimentally (see, e.g., Fig. 1 of Del Negro et al. 2002a). The model of fast excitatory synaptic dynamics used is as previously specified in Butera et al. (1999b). The results presented are from networks with all-to-all excitatory coupling. Similar results were obtained with sparsely connected networks (see following text).

Network simulations were run on Linux and Unix workstations using C code (Butera et al. 1999b) and the results were analyzed using Perl scripts and MATLAB. Network simulations were run for 2 min of simulation time with the first 30 s being ignored to allow start-up transients to decay.

Data analysis

An automated burst-detection technique was implemented by generating a combined histogram of spike times from every cell in the network (Fig. 4). The maximum and minimum amplitude of the histogram was calculated and the difference between those values was compared with a threshold. If the threshold was met and if the amplitude of the histogram remained <10% of the maximum amplitude for some minimum amount of time, then the output was defined as a burst. The values for amplitude threshold and minimum time <10% of maximum amplitude were chosen by visual inspection of the histogram data. This method classifies the results illustrated in Fig. 4, A and C as bursting, although the network producing the pattern illustrated in Fig. 4B is classified as a nonbursting network. For the results reported herein, only “regular” bursting patterns were considered; therefore the bursting patterns in Fig. 4, A and C were further analyzed. To determine whether regular bursting patterns were present, the burst period (BP), burst duration (BD), and burst amplitude were measured. Before measuring these values, the data were smoothed using a 20-point moving average. After smoothing, the start of a burst was calculated using a rising phase threshold of 30% of the maximum amplitude and the end of the burst was calculated using a falling phase threshold of 10% of the maximum amplitude. The BD is the time between these two thresholds, the BP is the time between the rising phase of two subsequent bursts, and the burst amplitude is the maximum amplitude measured during the BD. If each of these values had a coefficient of variation among all bursts in a trace <20% and if more than two bursts were found in the 90-s time window, then the bursting was defined as regular and is hereafter referred to as networkwide bursting. Thus the automated burst detection would define the bursts in Fig. 4C as irregular because the BP has a coefficient of variation >20%. The bursting pattern in Fig. 4A meets all requirements for regular networkwide bursting. This analysis scheme was verified in preliminary simulations studies by inspection of raster plots of the 50 cell-activity patterns to confirm the degree of cell synchrony between successive cycles. The values reported for the frequency of networkwide bursting were calculated using the inverse of the mean BP.

FIG. 4.

Network activity modes. A1C1: raster plots. A2C2: network activity, defined as histograms of spike times across the population (bin size = 10 ms). A1 and A2: regular bursting network. Most cells are bursting. B1 and B2: nonbursting network. Most cells are silent or tonically firing. C1 and C2: irregular bursting network. Burst period and amplitude are highly variable.

Once the existence of networkwide bursting has been established, input and output ranges are computed. Input range is defined as the size of parameter space where networkwide bursting occurs. The simulation data for input and output ranges are analyzed in some cases by grouping the number of PMs into bins where bin 0 is the special case when the network contains no PMs, bin 1 contains 1–5 PMs, bin 2 contains 6–10 PMs, and so forth. Input range is then defined as the percentage of the maximum number of simulations inside a bin that produce networkwide bursting, i.e., the fraction of the range of gtonic where networkwide bursting occurs. For example, bin 1 contains five values of PMs (from 1 to 5 PMs) and 16 values of gtonic (from 0 to 1.5 nS in 0.1-nS increments) for each number of PMs, giving 80 total simulations inside that bin. If networkwide bursting is found in only four of the 80 simulations (as is the case for bin 1 where gsyn = 0.2 nS; see following text), then the input range would be 4/80 or 5%. Output range is computed by subtracting the maximum frequency of networkwide bursting for a given number of PMs from the minimum frequency of networkwide bursting for the same number of PMs (as gtonic is varied) and averaging this range across all PMs in that bin. Therefore for a given number of PMs, if the network bursts with a maximum frequency of 0.75 Hz when gtonic = 0.9 nS and the network bursts with a minimum frequency of 0.25 Hz when gtonic = 0.1 nS, the output range is defined as 0.50 Hz for that number of PMs.



Figure 5 illustrates the existence and frequency of networkwide bursting as the number of PMs and excitatory drive (gtonic) are varied for five values of excitatory synaptic coupling (gsyn). The color indicates the frequency of networkwide bursting; black regions indicate the absence of networkwide bursting. The simulations in the black region could produce a network containing neurons that are all or mostly silent, a network containing neurons that are all or mostly tonically firing, or a network that produces an irregular bursting pattern. Each point in Fig. 5 is the result of a single simulation with a new randomly generated set of parameters. It is this independent parameter selection for each simulation that causes the “noise,” or variability among adjacent points. Figures 68 further quantify the results shown in Fig. 5; Figs. 9 and 10 illustrate the averaged results of multiple simulations using a single value for gsyn. Figure 11 provides a method of displaying our results that will allow a straightforward comparison to future experimental data.

FIG. 5.

Results of varying number of PMs (x-axis), gtonic (y-axis), and gsyn (from top left, 0.075 and 0.1 nS; from bottom left 0.15, 0.2, and 0.3 nS). Color indicates frequency of networkwide bursting. Black regions indicate the absence of networkwide bursting and could be either silence, a few cells bursting, irregular bursting, or tonic firing. Each point represents a distinct simulation with different randomly generated parameters.

FIG. 6.

Input range for several values of gsyn. Input range is computed as the percentage of the total number of simulations where networkwide rhythmic bursting occurs. x-axis is bin number where bin 0 is no PMs, bin 1 is 1–5 PMs, bin 2 is 6–10 PMs, and so forth.

Networkwide bursting

At lower levels of synaptic conductance (Fig. 5, top), more PMs are required to generate networkwide bursting. At higher levels of synaptic conductance (Fig. 5, bottom right), networkwide bursting is seen with fewer PMs, but the range of frequencies produced by the network is dramatically reduced. Other levels of synaptic conductance (>0.3 and <0.075 nS) were explored, but did not result in significant amounts of networkwide bursting. Also, increased excitatory drive (gtonic >1.5 nS) did not produce networkwide bursting at any level of synaptic conductance or for any fraction of PMs. At these high levels of excitatory drive, only tonically firing activity within the network was seen. The results of Fig. 5 also suggest that, for a fixed level of synaptic conductance and gtonic, reducing the number of PMs, in general, reduces the frequency of networkwide bursting until the rhythm is eventually abolished.

Input range and output range

Figure 6 is a plot of the input range for all values of gsyn as a function of K and levels of excitatory coupling. For the lower values of gsyn (0.075, 0.1, and 0.15 nS), the input range displays a generally increasing trend as the percentage of PMs increases from 0 to 100% of the population, with the lowest levels of synaptic conductance failing to produce any networkwide bursting until >20 or 25 PMs are in the population (bin 5, gsyn = 0.1 nS; bin 6, gsyn = 0.075 nS). For the larger values of gsyn (0.2 and 0.3 nS), the input range increases as the number of PMs is initially increased and then decreases slightly for higher percentages of PMs. Networks with gsyn values of 0.15 and 0.2 nS span the greatest range of input ranges (Fig. 6).

Table 2 quantifies the input and output ranges for three cases: a network containing no PMs, a network where less than (or equal to) half of the population is made up of PMs, and a network where more than half of the population is made up of PMs. Only two of the five synaptic conductances explored produced networkwide bursting when no PMs are present in the network (Fig. 6, bin zero). At gsyn values of 0.2 and 0.3 nS, networkwide bursting occurs with no PMs >13% of the input range. The average input range for a network with 1–25 PMs at gsyn values of 0.2 and 0.3 nS is 25 and 41%, respectively. The average input range for a network with more than half of its population made up of PMs at gsyn values of 0.2 and 0.3 nS is 52 and 39%, respectively (Table 2). Thus for a given level of synaptic conductance, increasing the number of PMs in the network can more than triple the input range compared with a network containing no PMs.

View this table:

Input and output ranges

Figure 7 illustrates the output range of the network for each level of synaptic conductance. The maximum and minimum frequencies for each number of PMs in a bin are given by solid and open circles, respectively. The frequencies measured in our simulations (from about 0.04 to 1.0 Hz) are similar to those measured by Del Negro et al. (2001; from about 0.05 to 0.8 Hz). For most cases, the output range increases as the number of PMs is increased (Fig. 7). Networks containing no PMs were able to burst (bottom center and right panels of Fig. 7); however, the output range in this case is considerably reduced compared with simulations with larger percentages of PMs. At gsyn values of 0.2 and 0.3 nS, the output range with no PMs is 0.20 and 0.04 Hz, respectively. The average output range when gsyn = 0.2 nS with 1–25 PMs and with 26–50 PMs is 0.19 and 0.43 Hz, respectively. When gsyn = 0.3 nS, the average output range with 1–25 PMs and with 26–50 PMs is 0.13 and 0.12 Hz, respectively (Table 2). The largest attainable output ranges for a given number of PMs with gsyn values of 0.15 and 0.2 nS are 0.88 Hz (with 46 PMs) and 0.71 Hz (with 41 PMs), respectively. These values are more than triple the largest output range obtained with no PMs in the network.

FIG. 7.

Output range for several values of gsyn. Results of varying number of PMs (x-axis is binned number of PMs), gtonic, and gsyn (from top left, 0.075, 0.1, 0.15, 0.2, and 0.3 nS). Solid line is frequency range, computed by subtracting the maximum frequency of networkwide bursting for a given number of PMs (closed circles) from the minimum (nonzero) frequency of networkwide bursting for the same number of PMs (open circles) and averaging this range across all PMs in that bin. x-axis is bin number where bin 0 is no PMs, bin 1 is 1–5 PMs, bin 2 is 6–10 PMs, and so forth, and y-axis is frequency.

Input–output range trade-offs

For a given level of gsyn, input range (Fig. 6) and output range (Fig. 7) were quantified and averaged across a range of K. Figure 8 plots the input range and the output range versus each other as K and gsyn are varied, allowing a compact visualization of how varying K and gsyn affect these metrics. Figure 8A is a plot of the input range versus output range for different numbers of PMs. The plot of Fig. 8A displays a clockwise trend as the level of synaptic conductance is increased for each number of PMs. In general, the input and output ranges increase as the number of PMs is increased. There is an optimal value of gsyn for a given number of PMs that maximizes the input and output ranges. Figure 8B is a plot of the input range versus output range for different levels of synaptic conductance. Here, the number of PMs increases as the trace progresses from left to right. Figure 8B emphasizes the depressing effect on output range of increasing synaptic strength. Weaker coupling produces a smaller input range and coupling also strongly produces a smaller output range. The moderate levels of synaptic conductance that provide the largest input and output ranges are found with gsyn values of 0.15 and 0.2 nS. At these moderate levels of synaptic conductance, the largest gain (or largest rate of change) in input and output ranges occurs when the number of PMs increases to >20–40% of the population (see, for example, the large increase between the second and third points where gsyn = 0.2 nS in Fig. 8B; see also Fig. 5).

FIG. 8.

Input range vs. output range for different numbers of PMs (A) and for 5 different levels of synaptic coupling (B). A: synaptic conductance increases from 0.075 to 0.3 nS as each trace progresses clockwise. B: number of PMs increases as each trace progresses from left to right. Numbers of PMs are grouped according to the legend of A (i.e., 0 PMs, 1–10 PMs, etc.).

Multiple simulations

Because the optimal synaptic conductance level was found to be between gsyn = 0.15 nS and gsyn = 0.2 nS, four additional simulations were run at both levels of synaptic conductance. This was to validate general trends just described and average over the variability of individually randomly generated simulations. The input and output ranges for each simulation are given in Fig. 9 and the average of all five simulations is plotted as a bold line. The input range is fairly consistent among all simulations for gsyn values of 0.15 and 0.2 nS (Fig. 9, A and B). The output range for gsyn = 0.2 nS is more variable (Fig. 9D) but both display a clear upward trend as the number of PMs is increased (Fig. 9, C and D). Previously, the largest output range obtained for a network with no PMs was 0.20 Hz when gsyn = 0.2 nS. However, none of the four additional simulations at gsyn = 0.2 nS produced networkwide bursting at more than a single value of gtonic with no PMs in the network. This reduces the average output range with no PMs in the network to 0.04 Hz.

FIG. 9.

Input (A and B) and output (C and D) range of 5 simulations where gsyn = 0.15 nS (A, C) and 0.2 nS (B, D). Bold line is the average of 5 runs.

Figure 10, A and B displays the average results of varying gtonic and number of PMs for all simulations with gsyn values of 0.15 and 0.2 nS. When computing the average, simulations that did not produce networkwide bursting were ignored (rather than including a zero frequency in this average). Averaging multiple runs removes much of the variability among adjacent points seen in the results from a single run (compare Fig. 10, A and B to bottom left and center panels of Fig. 5). Figure 10, C and D reveals the likelihood of choosing a random set of parameters that can produce networkwide bursting with a certain number of PMs in the population for a given level of excitatory drive. In Fig. 10, C and D, the (black and white) color bar indicates in how many of the five simulations networkwide bursting was measured. The white regions are where bursting occurs for all five simulations and the black regions are where no networkwide bursting is recorded. Figure 10, B and D illustrates that the network is able to burst with fewer numbers of PMs and is also capable of producing a large range of frequencies with fewer PMs. However, the likelihood of producing networkwide bursting in this region is minimal (Fig. 10D, top left). The network does produce bursting when few PMs are present, but it does not consistently produce bursting at the same points in parameter space. This is the cause of the variability seen in Fig. 9D.

FIG. 10.

A and B: average results of varying number of PMs (x-axis) and gtonic (y-axis) for 5 simulations where values of gsyn are 0.15 nS (A) and 0.2 nS (B). Individual simulations where no bursting was found were not included in the average. C and D: likelihood of bursting in 5 simulations where values of gsyn are 0.15 nS (C) and 0.2 nS (D). White region indicates bursting for all 5 simulations and black region indicates no bursting for any simulation.

Sparsely connected networks

The simulations were repeated using a sparsely connected population and the results were compared with the all-to-all coupling simulations given earlier. For these simulations, the probability of connection was reduced to 10% (as an example of very sparse connectivity) and the level of synaptic conductance was increased so that the mean level of synaptic current provided to each cell remained roughly the same. For the sparsely connected networks, we explored three values of synaptic conductance, gsyn values of 1.0, 1.5, and 2.0 nS, that roughly correspond to gsyn values of 0.1, 0.15, and 0.2 nS using the all-to-all connectivity presented earlier. The results of the simulations using the sparsely connected networks (not shown) are completely consistent with all results using all-to-all connectivity. Both input and output ranges of the sparsely connected networks increased as the number of PMs increased in a similar fashion to the corresponding levels of synaptic conductance with all-to-all coupling. Also, the magnitudes of the ranges of frequencies measured were similar.


These simulations were performed to investigate the functional significance of PM cells in the pBC excitatory network. The simulations use two types of heterogeneity: cell-type heterogeneity (PMs and NPMs) and parameter heterogeneity within a given cell type that is based on experimental data. Before running the data-based simulations presented here, we performed several simulations (not shown) using less heterogeneous parameter distributions: a semihomogeneous parameter space (gLeak was kept constant for all cells; gNaP was constant within each subpopulation of PMs or NPMs) and a one-dimensional heterogeneous parameter space (gLeak was kept constant for all cells whether PM or NPM; mean gNaP was different for PMs and NPMs and varied with SD of ±10%, and then ±30%). The semihomogeneous simulations produced networks with reduced input and output ranges compared with networks with more heterogeneous parameter distributions, which might be expected based on previous studies examining the role of heterogeneity in this network (Butera et al. 1999b). Although there were minor quantitative differences, both the simulations using less heterogeneous and less data-based parameter distributions and the simulations using sparsely connected networks provide results that are qualitatively similar to the results of the data-based simulations presented here. The major conclusion of this study was reiterated in every simulation performed, that is, pacemakers can increase the input and output ranges of a bursting network.


The absolute values for input range are not meaningful because the percentages calculated are a function of the size of parameter space chosen (gtonic varies from 0 to 1.5 nS in 0.1-nS increments). If gtonic were allowed to increase to ≤2.0 nS, then all calculated percentages would be reduced. However, the relative differences between the input ranges for different numbers of PMs and for different levels of gsyn are important. The input ranges given in Fig. 6 and Table 2 reveal that networkwide bursting is easier to obtain when the fraction of the network population made up of PMs is increased. Like input range, output range shows a generally increasing trend as the number of PMs is increased. This highlights an important difference between the controllability potential of synaptic conductance versus the number of PMs in a network. Increasing gsyn can increase input range but at the cost of output range. Alternatively, increasing the number of PMs in the network can increase both input range and output range.

Comparing both the input range and the output range for several values of gsyn reveals an optimal synaptic conductance near gsyn values of 0.15 and 0.2 nS. Therefore we ran multiple simulations at these levels of synaptic conductance. Examining the results of these simulations revealed that large output ranges were possible with few PMs at gsyn = 0.2 nS, but the output range in this region was highly variable because of the inability to unfailingly generate networks capable of producing networkwide bursting. In general, the network input–output range increased as the percentage of PMs increased under these optimal synaptic coupling conditions.

The reported estimates for the number of PMs in the pBC vary from 5 to roughly 50% for different in vitro slice preparations and recording conditions (Del Negro et al. 2002a, 2005; Koshiya and Smith 1999; Pagliardini et al. 2005; Pena et al. 2004). Because the PM or NPM state of a pBC cell is modulatable (Pena and Ramirez 2002, 2004), it is possible that the number of PMs in this network is modulated to meet particular dynamic range requirements. Indeed, our current and original model analyses indicate that any neuromodulatory conditions that change the gNaP/gLeak ratios within the network can potentially change the relative numbers of PMs and NPMs and thus the network's dynamic performance. Thus the varying estimates for the number of PMs in the pBC under different recording conditions may in part reflect the modulated state of the pBC. With optimal values of gsyn in our simulations, the largest increase in dynamic range occurs when the number of PMs increases to ≳20% of the population (see Figs. 810). We are currently using reduced-order models to investigate the mechanism and significance of this threshold-like effect of the number of PMs on network dynamics.

Limitations of the model

These simulations use a minimal model of pBC cells that lack some voltage-gated ionic currents known to exist in these cells, such as the low-voltage–activated or high-voltage–activated calcium channels (Elsen and Ramirez 2005). This minimality undoubtedly plays a role in the inability of the model to perfectly match the experimental data. In constructing our network models, we sought to rigorously incorporate conductance density distributions, their variability, and voltage dependencies that are consistent with the experimental data. An important model parameter that has not yet been directly verified by measurements, however, is the kinetics for the voltage-dependent slow inactivation of NaP incorporated in the model. This parameter controls the dynamics of the single-neuron rhythmic bursting cycle and synchronized rhythmic bursting at the population level.

Because model parameters were chosen from a two-dimensional heterogeneous parameter space, ideally, multiple simulations should be run so that any variability among simulations can be averaged out. This is computationally time consuming. Figure 5 contains almost 4,000 simulations, where each simulation is a network of 50 cells. However, our results indicate that more simulations are not necessary. When multiple simulations were run for a single value of synaptic conductance, all simulations displayed similar trends. Also, as previously stated, all simulations using less heterogeneous parameter distributions and sparsely connected networks provided results consistent with the results presented in this study.

The simulations that produced networkwide bursting with few or no PMs typically provided minimal output ranges. Although some simulations with few PMs were able to produce large output ranges (Fig. 9D), the likelihood of producing networkwide bursting here was low (Fig. 10D). Attempts were made to choose data-based parameters for the simulations. However, the limitations of the minimal model and of the parameter selection technique could be the source of the inability to unfailingly produce robust networkwide bursting with fewer numbers of PMs in the population.

We should note that, although PMs are not required to produce networkwide bursting in our model (Butera et al. 1999b), the NaP current is required. If this slowly inactivating inward current were completely removed from our model, there would be no slow process in the network to terminate networkwide bursting activity once it has begun. Thus we emphasize that NaP provides not only a mechanism for rhythmic bursting of individual PMs (when isolated from synaptic input), but also a mechanism for synchronized rhythmic bursting and its termination at the population level. This mechanism results from the voltage dependency and kinetics of NaP inactivation and the dynamic interaction of NaP current with phasic excitatory synaptic drive currents in the coupled population of cells (see Fig. 10 of Butera et al. 1999b).

Whether the NaP current provides the primary mechanism for burst termination and regenerative reactivation of the network in in vitro slice preparations remains to be clearly resolved. Published investigations of this mechanism appear to appreciably depend on the experimental preparation (i.e., thin vs. thick slices in vitro vs. the more intact network in situ) (Del Negro et al. 2005; Paton et al. 2006; Pena et al. 2004; Rybak et al. 2003; Tryba et al. 2006). Our results (Fig. 5) suggest that at least for preparations with a relatively isolated pBC network (e.g., thin in vitro slice preparations) and for a fixed level of synaptic conductance, reducing gNaP and thus the number of PMs, in general, reduces the frequency of networkwide bursting until the rhythm is eventually abolished. This agrees with some experimental results where gNaP was reduced pharmacologically by riluzole or low concentrations of tetrodotoxin (TTX), presumably reducing gNaP/gLeak ratios to parameter space regions with few PMs (Koizumi and Smith 2004; Rybak et al. 2003). Our model does not include the cadmium-sensitive bursters found in the pBC of slice preparations from mice (Del Negro et al. 2005; Pena et al. 2004; Thoby-Brisson and Ramirez 2001; Tryba et al. 2006), but not rats (Del Negro et al. 2005). Further modeling and simulations are required to determine the effects of adding these pacemakers to our network models. At present, little is known about the biophysical parameters of the CaN current and these cells may (Thoby-Brisson and Ramirez 2001) or may not (Del Negro et al. 2005) provide voltage-dependent control of PM oscillation frequency, which would impart a very different contribution to the input–output range of the network. Indeed, a unique property of NaP-dependent mechanisms may be the enabling of a broad range of tunable network oscillation frequencies, which is inherent in the voltage-dependent properties of this conductance mechanism. In considering the functional significance of any type of PM mechanism, the issue of frequency control may be fundamental. That is, it is not simply the capability of a particular cellular pacemaker mechanism when coupled with excitatory synaptic interactions to organize a networkwide rhythm in the pBC, but it is the inherent ability of any mechanism to allow functionally for frequency control over a wide dynamic range. This control is a prominent feature of experimental data, at least for the isolated pBC in slices in vitro under conditions where tonic drive is varied (e.g., Del Negro et al. 2001).

Rhythmic networks

Although we have analyzed dynamics of excitatory networks incorporating bursting cells to investigate rhythm generation mechanisms in the pBC, our results would pertain to rhythmogenesis in the isolated pBC as in slice preparations in vitro, or under hypoxic conditions in perfused brain stem preparations in situ (Paton et al. 2006), where there is evidence for NaP-dependent mechanisms. As we originally pointed out (Smith et al. 2000), when the pBC is embedded in the brain stem respiratory network in more intact states of the system, as in vivo, phasic (and possibly tonic) inhibitory inputs to the pBC become important for dynamically regulating the evolution/termination of pBC inspiratory activity (Tryba et al. 2003). We proposed (Smith et al. 2000), for example, that network-based phasic inhibition contributes importantly to the termination of pBC networkwide bursting in the intact system, as opposed to termination solely by the slow inactivation kinetics of NaP and interactions of NaP and leak currents. Thus it remains an important theoretical and experimental problem to understand rhythm generation when inhibitory network mechanisms are overlaid on the excitatory network interactions that we have analyzed here.

Although our model is based on intrinsic bursting cells in the pBC, this study is representative of a particular class of rhythmic bursting networks—those where the connectivity is predominantly mediated by excitatory synapses. Examples include rhythmic bursting in transverse brain stem slices (Koshiya and Smith 1999; Smith et al. 1991), rhythmic activity in disinhibited embryonic spinal cords (O'Donovan 1989; Streit 1993), rhythmic activity in disinhibited medial septum and diagonal band complex (Manseau et al. 2005), and 7- to 14-Hz oscillations in motor cortex (Castro-Alamancos et al. 2006). In such networks, the synaptic excitation is responsible for the coordination and spread of activity through the network at the start of an episode of activity; in other cases such excitation may play a key role in the initiation of a burst of activity as well. In all these cases, burst episode initiation is attributed to some combination of intrinsic cellular properties and/or recurrent excitatory coupling. Burst termination can be explained by the slower kinetics of ionic currents intrinsic to the component neurons (e.g., this manuscript) or the slower kinetics associated with synaptic depression (Tabak et al. 2000).

Model predictions

To date, we are unaware of any experimental studies that quantify dynamic range under conditions that putatively alter the number of PMs. Based on the simulation results, we expect both the input range and the output range of the network would be reduced as the fraction of PMs in the network is reduced. Likewise, if the conductance density of NaP channels is augmented, then we predict that neurons that were previously NPMs would become PMs and would produce an increased dynamic range. Some have claimed (Del Negro et al. 2002b, 2005) and our modeling studies have shown that networkwide bursting is possible after sufficiently reducing NaP such that no PM activity is seen in the network (Butera et al. 1999b). When PMs are blocked experimentally, an increased level of depolarization is required to restore rhythmic activity (Del Negro et al. 2005). Our simulations confirm these results. If no (or few) PMs are present in the network, increased levels of gtonic are required to produce networkwide bursting (Figs. 5 and 10B). However, our simulations suggest that under these conditions the dynamic range of the network is greatly reduced. Substance P was previously shown to reactivate networkwide bursting after ostensibly blocking all PMs (Del Negro et al. 2005; Tryba et al. 2006). Even if this reactivation is caused by the apparent ability of substance P to temporarily restore PM activity to some cells where PM activity was previously abolished (Tryba et al. 2006), based on our simulation results, we would still expect the dynamic range of the network to be reduced under these conditions. Experiments to test the model predictions would involve measuring the range of frequencies produced by the network while varying, for example, the level of tonic activation of α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) receptors or tonic excitation with extracellular potassium concentration under control conditions and when gNaP and thus PM activity have been attenuated by application of riluzole or small amounts of TTX. Our model analog to the experimental figure this would produce is given in Fig. 11, where we plot frequency versus gtonic for different numbers of PMs. Experimentally varying the level of AMPA receptor activation or extracellular potassium within the pBC is qualitatively equivalent to varying EK (Del Negro et al. 2001; Rybak et al. 2003), EL (Butera et al. 1999a), or gtonic (Butera et al. 1999b) in the model. Figure 11 can be compared with future experimental tests of these model predictions. In Fig. 11, the nonmonotonic increase in frequency as gtonic is increased is caused by heterogeneity among simulations. This variability is reduced when averaging the results of multiple simulations (see bold lines in Fig. 9).

FIG. 11.

Frequency vs. gtonic for different numbers of PMs when gsyn = 0.15 nS. Reducing the number of PMs reduces both the input range and the output range of the network. No bursting is produced for 0 PMs or 1–10 PMs.

Hypoxia increases NaP conductance density (Hammarstrom and Gage 1998, 2002; Horn and Waldrop 2000) and recent studies showed the importance of NaP-dependent bursting cells during hypoxia (Paton et al. 2006; Tryba et al. 2006). If hypoxia increases the number of PMs in the network, this would cause a shift in parameter space to a region capable of producing higher-frequency rhythms. This is a potential mechanism for obtaining the initial increase in burst frequency seen during hypoxia (Solomon et al. 2000; Telgkamp and Ramirez 1999; Tryba et al. 2006). Again, the dynamic range of the network could be examined at different levels of pBC hypoxia. In general such comparisons under any experimental conditions that are found to either selectively augment gNaP or augment gNaP/gLeak ratios to increase the number of PMs would be instructive.

In conclusion, although the significance and abundance of PMs in the pBC are not completely understood, the existence of these cells is not questioned. We used modeling studies to quantitatively explore the role that one of the main types of PMs (i.e., NaP-dependent) found in the pBC plays in this bursting network. Networks containing no PMs were able to produce regular, synchronous networkwide bursting (albeit with a low probability of bursting, as well as low input and output ranges), demonstrating that PMs are not critical for rhythm generation, provided that excitatory synaptic coupling strength is high. Thus the NaP current can provide not only a mechanism for rhythmic bursting of individual PMs, but also a mechanism for synchronized rhythmic bursting and its termination at the network level. However, our modeling studies suggest that including PMs in the network allows the input and output ranges to more than triple compared with networks with no PMs. Indeed, the fraction of PMs profoundly affects the controllability of the rhythm. Unlike synaptic coupling that requires a trade-off between input range and output range, increasing the number of PMs can increase both the input range and the output range of the network. Additional experimental tests as suggested by our results will be required to either confirm or refute our model prediction that PMs increase the robustness of rhythm generation.


The model used for all simulations is a modified version of Model 1 of Butera et al. (1999a) and is based on a single-compartment Hodgkin–Huxley (HH) formalism. The membrane potential is found using the differential equation Math where V is the membrane potential (mV), Cm is the whole cell capacitance (21 pF), t is time (ms), Itonic is the excitatory drive current (nA), Isyn is the synaptic current (nA) from other pBC cells in the network, and Iionic are the ionic currents listed below and have the form Math where Eionic is the equilibrium reversal potential for the ionic species carried by the current and Math where is the maximum conductance of each current, and x is the product of one or more gating variables raised to integer powers, as subsequently described.

The dynamics of the conductances of the ionic currents regulated by voltage-dependent activation or inactivation variables are described according to Math Math Math where x(V) is the steady-state voltage-dependent (in)activation function of x and τx(V) is the voltage-dependent time constant. x(V) is a sigmoid with a half-(in)activation at V = θx and a slope factor σx. τx(V) is a bell-shaped curve that has a maximal value τ̄x at V = θx and a half-width determined by σx. Thus each gating variable is described by only three parameters.

Action potentials in the model are generated by a fast Na+ current (INa) and a delayed-rectifier K+ current (IK). The equations for these currents are inspired by the HH formulation, but the gating variables satisfy a reduced voltage-dependent description Math Math where the parameters of INa are Na = 28 nS, ENa = 50 mV, θm = −34 mV, and σm = −5 mV; the parameters of IK are K = 11.2 nS, EK = −85 mV, θn = −29 mV, σn = −4 mV, and τ̄n = 10 ms.

The persistent sodium current is described by the following equation Math where ENa = 50 mV, θm = −45.1 mV, σm = −5 mV, θh = −53 mV, σh = 6 mV, and τ̄h = 10,000 ms. NaP is varied in each simulation as described earlier.

The model has a passive leak current defined as Math This current is K+ dominated. EL is set to −70 mV and Leak is varied in each simulation as previously described. The reader is referred to Butera et al. (1999a,b) for a more detailed description of the model.


This work was supported by National Institute of Mental Health Grant R01-MH-62057 and, in part, by the Intramural Research Program of the National Institute of Neurological Disorders and Stroke, National Institutes of Health.


  • 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