Models of Respiratory Rhythm Generation in the Pre-Bötzinger Complex. III. Experimental Tests of Model Predictions

Christopher A. Del Negro, Sheree M. Johnson, Robert J. Butera, Jeffrey C. Smith


We used the testable predictions of mathematical models proposed by Butera et al. to evaluate cellular, synaptic, and population-level components of the hypothesis that respiratory rhythm in mammals is generated in vitro in the pre-Bötzinger complex (pre-BötC) by a heterogeneous population of pacemaker neurons coupled by fast excitatory synapses. We prepared thin brain stem slices from neonatal rats that capture the pre-BötC and maintain inspiratory-related motor activity in vitro. We recorded pacemaker neurons extracellularly and found: intrinsic bursting behavior that did not depend on Ca2+ currents and persisted after blocking synaptic transmission; multistate behavior with transitions from quiescence to bursting and tonic spiking states as cellular excitability was increased via extracellular K+concentration ([K+]o); a monotonic increase in burst frequency and decrease in burst duration with increasing [K+]o; heterogeneity among different cells sampled; and an increase in inspiratory burst duration and decrease in burst frequency by excitatory synaptic coupling in the respiratory network. These data affirm the basis for the network model, which is composed of heterogeneous pacemaker cells having a voltage-dependent burst-generating mechanism dominated by persistent Na+ current (I NaP) and excitatory synaptic coupling that synchronizes cell activity. We investigated population-level activity in the pre-BötC using local “macropatch” recordings and confirmed these model predictions: pre-BötC activity preceded respiratory-related motor output by 100–400 ms, consistent with a heterogeneous pacemaker-cell population generating inspiratory rhythm in the pre-BötC; pre-BötC population burst amplitude decreased monotonically with increasing [K+]o (while frequency increased), which can be attributed to pacemaker cell properties; and burst amplitude fluctuated from cycle to cycle after decreasing bilateral synaptic coupling surgically as predicted from stability analyses of the model. We conclude that the pacemaker cell and network models explain features of inspiratory rhythm generation in vitro.


Rhythmic breathing movements in mammals are hypothesized to originate from patterns of neural activity generated in the pre-Bötzinger complex (pre-BötC), a specialized region of the ventrolateral medulla (Gray et al. 1999; Smith et al. 1991). Neurons in the pre-BötC are both necessary and sufficient to generate inspiratory-related motor output in vitro (Rekling and Feldman 1998; Smith et al. 1991), and perturbations or lesions of this region disrupt inspiratory activity in vivo (Hsieh et al. 1998; Koshiya and Guyenet 1996; Ramirez et al. 1998; Solomon et al. 1999). Therefore the oscillatory network contained in the pre-BötC putatively represents the most rudimentary substrate or kernel for generation and regulation of respiratory rhythm (at least in vitro) and can be retained in thin slice preparations from neonatal rodents that generate inspiratory-related motor activity.

Respiratory rhythm generation does not require synaptic inhibition (Feldman and Smith 1989; Gray et al. 1999), and a subset of inspiratory interneurons in the pre-BötC are bursting-pacemaker neurons synchronized by non-N-methyl-d-aspartate (NMDA) fast excitatory synapses (Johnson et al. 1994; Koshiya and Smith 1999a; Smith et al. 1991; Thoby-Brisson and Ramirez 2000; Thoby-Brisson et al. 2000). These data suggest that the rhythm-generating mechanism in vitro incorporates an excitatory network of synaptically coupled pacemaker neurons (for review, see Rekling and Feldman 1998).

In the first two papers of this series, Butera et al. created mathematical models of inspiratory pacemaker neurons (Butera et al. 1999a), which were assembled to form a network model of the rhythm-generating kernel (Butera et al. 1999b). These models posited burst-generating mechanisms for pacemaker neurons and examined how cellular heterogeneity, excitatory synaptic coupling, and tonic excitation could influence network-level rhythm generation as well as the behavior of individual cells in the context of network activity. We have now evaluated the basis for the models and the testable predictions that emerged from the modeling studies at cellular and network levels. We found that the pacemaker neuron and network models successfully predicted many behaviors observed in vitro. These new data affirm many of Butera et al. (1999a,b)'s theoretical conclusions regarding pacemaker cell behaviors and the roles of excitatory synapses and cellular heterogeneity in generation and control of respiratory rhythm in vitro. Several deficiencies of the models are also identified here and we suggest extensions of the models to refine our understanding of rhythm generation.


Experimental methods

Thin transverse slices (350 μm-thick) with rostral and caudal ends bordering the pre-BötC were cut from the medulla of neonatal rats (P0–P3) in artificial cerebrospinal fluid (ACSF) containing (in mM) 128.0 NaCl, 3.0 KCl, 1.5 CaCl2, 1.0 MgSO4, 21.0 NaHCO3, 0.5 NaH2PO4, and 30.0d-glucose, equilibrated with 95% O2-5% CO2 (27°C, pH = 7.4), as originally described (Smith et al. 1991). Slices were pinned down in a 2-ml recording chamber and perfused with ACSF at 5 ml/min. Rhythmic respiratory activity was maintained by raising the ACSF K+ concentration ([K+]o) to 5–8 mM. Inspiratory-related motor discharge (Smith et al. 1990) was recorded from the hypoglossal nerve (XIIn) rootlets (also captured in the slice) using fire-polished glass suction electrodes (60–90 μm ID) and a Cyberamp 360 (Axon Instruments, with variable gain and a 0.3- to 1-kHz band-pass filter (Fig.1). In many experiments, inspiratory neuron population activity was simultaneously recorded locally in the pre-BötC using “macropatch” suction electrodes (∼100 μm ID; Fig. 1). The slices were typically cut so that the caudal end of the pre-BötC was exposed, and the population recordings were made from this caudal surface. Both the XIIn and pre-BötC recordings were rectified and smoothed by analog integration and acquired digitally with raw signals at 4 kHz in PowerLab (ADInstruments, Inspiratory neurons were recorded extracellularly in the pre-BötC using 0.5-M sodium acetate-filled microelectrodes (8–12 MΩ).

Fig. 1.

A schematic showing the 350-μm-thick neonatal rat medullary slice preparation and typical recording configuration. Population-level recordings obtained using “macropatch” suction electrodes applied locally in the pre-Bötzinger complex (pre-BötC) and from roots of the XIIn were acquired and displayed as raw signals and after rectification/smoothing with a leaky analog integrator (∫pre-BötC and ∫XII). Sample traces show consecutive inspiratory bursts recorded from the pre-BötC and XIIn at 8-mM [K+]o and cycle-triggered averages of burst activity (average of 60 consecutive respiratory cycles) on an expanded time scale.

To identify pacemaker neurons, the excitatory synaptic transmission critical for respiratory network function in vitro (Funk et al. 1993; Ge and Feldman 1998; Koshiya and Smith 1999a) was blocked using 20-μM 6-cyano-7-nitroquinoxaline-2,3-dione disodium (CNQX). All sources of chemical synaptic transmission were blocked using low-Ca2+ ACSF containing (in mM) 120 NaCl, 8–12 mM KCl, 0.2 CaCl2, 1 MgSO4, 4 MgCl2, 21 NaHCO3, 0.5 NaHPO4, and 30 mM d-glucose (27°C, pH = 7.4).

In some experiments, the slice preparation was severed along the midline resulting in two bilaterally separated symmetrical “split slices” (Fig. 10 A); inspiratory activity was monitored simultaneously bilaterally via local recordings in the pre-BötC.

Experimental data analysis

Measurements of cellular and network activity such as inspiratory burst period, frequency, and duration were determined off-line using automated algorithms, hand-checked for accuracy. Inspiratory bursts obtained from XIIn or local pre-BötC recordings were detected by threshold crossings. We constructed baseline noise histograms (105 points) from quiescent intervals between bursts (i.e., the network expiratory phase) and fit Gaussian functions to the baseline noise distribution to determine the standard deviation. The event threshold was set at 2 SD greater than mean baseline noise, which ensures that the probability of selecting spurious events was P < 0.05. At high levels of excitability where net inspiratory discharge declined (see Fig. 9), we used fast Fourier transforms to corroborate the mean frequency calculations from the time domain. The threshold-crossing criteria applied to rhythmic XIIn discharge were used to cycle-trigger running averages of XIIn and pre-BötC activity (Figs. 8 and 9).

The action potentials and bursts from pacemaker neurons greatly exceeded baseline noise (e.g., Figs. 4 and 5) and were also selected by threshold-crossing algorithms. The first spike of four that occurred within a sliding 133-ms window (30 Hz) defined burst onset. The last action potential in the burst was determined from interspike intervals (ISIs) if two criteria were satisfied: the last spike preceded an ISI ≥500 ms and the next spike marked subsequent burst onset (defined by the onset criteria). Action potentials that occurred in the sliding window but failed to satisfy onset criteria were ignored. These criteria were also applied to model data and allowed us to distinguish inspiratory bursts from low-frequency spiking activity that sometimes occurred between inspiratory cycles. Burst duration was defined as the time spanning burst onset to offset.

All data presented as burst frequency in this study were first analyzed by computing burst period, defined as the interval from onset to onset in two consecutive bursts. Statistical tests were performed using periods before plotting as frequency (reciprocal of period) to avoid error that can arise from prior conversion to frequency. The changes in burst period and burst duration of single cells before and after CNQX application and/or low-Ca2+ solutions were assessed using paired t-tests.

Mathematical modeling and numerical methods


Butera et al. (1999a) modeled inspiratory pacemaker neurons of the pre-BötC using the fewest number of state variables and parameters needed to reproduce intracellular data. In this study, we used their model 1 (see Butera et al. 1999a for discussion of why model 1 is favored), where a fast-activating, slowly inactivating persistent Na+ current (I NaP) primarily constitutes the burst-generating mechanism. Other currents include Hodgkin-Huxley-like Na+ current (I Na) and delayed-rectifier-like K+ current (I K) for action potential generation, K+- and Na+-dominated leakage currents (I L(K) andI L(Na)), and tonic excitatory synaptic current (I tonic(e)). In published classification schemes for modeled bursting neurons, this model conforms to type I (Bertram et al. 1995) or fold-homoclinic (Izhikevich 2000).

Model 1 has three state variables. V Mand n are fast changing, whereV M is membrane potential (in mV) andn is a gating variable that activatesI K. The third variable h is a gating variable for slow inactivation ofI NaP and influences phase transitions during bursting. State variables evolve according to nonlinear ordinary differential equationsCV˙M=INaPINaIKIL(K)IL(Na)Itonic(e) Equation 1 n˙=(n(VM)n)/τn(VM) Equation 2 h˙=(h(VM)h)/τh(VM) Equation 3 Equation 1 describes membrane potential trajectory, where C is whole cell capacitance (in pF). Gating variablesn and h converge to steady statesn (V M) andh (V M) with kinetics determined by τn(V M) and τh(V M). The voltage-dependent functions for (in)activation and time constants arex(VM)={1+exp[(VMθx)/ςx]}1,x=n,h Equation 4 τx(VM)=τ¯x/cosh[(VMθx)/(2ςx)],x=n,h Equation 5wherex (V M) are sigmoidal activation functions with slopes proportional to 1/ςx at V M = θx and τx(V M) are bell-shaped functions with peak τ¯ xat V M = θxand width ∝ 2ςx. Action potential currentsI Na andI K are described byINa=g¯Nam3(VM)(1n)(VMENa) Equation 6 IK=g¯Kn4(VMEK) Equation 7where g¯ Na and g¯ K are maximal conductances (in nS),m (V M) is voltage-dependent activation, andE Na andE K are reversal potentials.I Na activation occurs on a much faster time scale than changes in membrane potential (i.e., M) due in part to the scaling of voltage change by membrane capacitance. ThereforeI Na is considered to activate instantaneously (Butera et al. 1999a). Variablen simultaneously activatesI K and inactivatesI Na via the (1 − n) term since these processes have similar voltage dependence and kinetics (Butera et al. 1999a).I NaP is described byINaP=g¯NaPp(VM)h(VMENa) Equation 8with maximal conductance g¯ NaPand instantaneous activationp (V M).

In this study, we replaced the original leakage current (I L) with specific K+ and Na+ components,I L(K) andI L(Na), described byIL(K)=gL(K)(VMEK) Equation 9 IL(Na)=gL(Na)(VMENa) Equation 10where g L(K) andg L(Na) are nongated conductances. This change enabled us to compare our experimental perturbation (changing [K+]o) toE K in the model; both depolarize neurons by modifying the driving force for K+.I tonic(e) models tonic excitatory inputItonic(e)=gtonic(e)(VMEsyn(e)) Equation 11where g tonic(e) is synaptic conductance and E syn(e) is the reversal potential of non-NDMA excitatory amino-acid (EAA) receptors. Single-cell simulations used g tonic(e)= 0, but for network simulationsg tonic-e ≠ 0 (Table1).

View this table:
Table 1.

Parameter distributions used for pacemaker-network simulations

Standard parameters for model 1 are: C = 21 pF, g¯ Na = 28 nS, g¯ K = 11.2 nS, g¯ NaP = 2.4 nS,g L(K) = 2.4 nS,g L(Na) = 0.4 nS,E Na = 50 mV,E K = −85 mV,E syn(e) = 0 mV, θm = −34 mV, ςm = −5 mV, θn = −29 mV, ςn= −4 mV, θp = −40 mV, ςp = −6 mV, θh = −48 mV, ςh = 6 mV, τ¯ n = 10 ms, τ¯ h = 104 ms, and g tonic(e) = 0 nS.


The respiratory rhythm-generating kernel was modeled as Nheterogeneous pacemaker neurons, coupled by non-NMDA fast excitatory synapses (Butera et al. 1999b). Phasic synaptic current (I syn(e)) was incorporated intoEq. 1 for cells of the network.I syn(e) in neuron j is the sum of excitatory synaptic input from N − 1 non-j cells in the population (all-to-all coupling)Isyn(e),j=i=1N1g¯syn(e),i,jsi(VMjEsyn(e)) Equation 12 s˙i=[(1si)s(VMi)krsi]/τs Equation 13where g¯ syn(e),i,jis synaptic conductance between neurons i and j, and si is the synaptic gating variable. Presynaptic action potentials activatesi, whose kinetics approximate fast excitatory synapses (τs = 5 ms andk r = 1) in respiratory neurons (Funk et al. 1993; Ge and Feldman 1998). The functions (V M) determines the steady-state postsynaptic receptor activation based on presynaptic membrane potential in neuron i. All simulations use N = 50 (see Butera et al. 1999b for choice of N).

To incorporate heterogeneity in pacemaker neurons the parameters g¯ NaP, g¯ L(K), and g¯ syn(e) were randomly assigned from normal distributions (Table 1). Sometimes (Figs. 7-9) we included afollower cell population of interneurons one excitatory synapse downstream from the rhythm-generating network. Follower cells do not synaptically project to cells in the kernel, nor to each other, and do not participate in rhythm generation. Follower cells were constructed using model 1 with g¯ NaP = 0 nS (they are nonpacemakers) and different coupling strength (Table1). The probability of synaptic connection between pacemaker neurons and follower neurons was 0.5. To simulate split-slice conditions, we divided the model into two independent (synaptically uncoupled) halves each containing N/2 cells. All other parameters in split-slice simulations remained the same.

Numerical methods

Computer simulations of the isolated pacemaker model used the CVODE numerical integrator (Scott D. Cohen and Alan C. Hindmarsh, and XPP software (Bard Ermentrout, Network simulations were coded in the C programming language and run on Pentium-Linux (Dell, and Ultrasparc-Solaris (Sun Microsystems, workstations, using a fifth-order Runge-Kutta-Fehlberg integration method with Cash-Karp parameters and adaptive time step (initial conditions randomly assigned) (Butera et al. 1999b;Press et al. 1992). Before collecting model data, 90 s of simulated time was allowed for initial transient decay. Network activity was displayed as a running histogram (adjustable bin size: 10–100 ms) of spike times computed across the network or as a raster plot of spike times in the population (Figs. 7 and 8).


Evaluation of pacemaker neuron model


Butera et al. (1999a) proffered model 1 as a minimal mathematical model of voltage-dependent bursting pacemaker neurons in the pre-BötC. Bursting is influenced by the excitability parameters that control baseline membrane potential (V M) such as applied current (I app) orE K, which acts via the K+-dominated leakage current. Here we selectedE K for the adjustable excitability parameter to more accurately compare our simulations with in vitro experiments where [K+]owas used to control excitability (see methods). At hyperpolarized V M, the model is quiescent (e.g., V M = −56 mV atE K = −80 mV, Fig.2 A). Bursting oscillations, the periodic alternation between bursts of action potentials and quiescent interburst intervals, emerge asE K is used to depolarize baselineV M. At the highestE K, the neuron exhibits tonic spiking (e.g., E K = −72 mV, Fig.2 A). Therefore the neuron is a “conditional” pacemaker since its oscillatory activity depends on voltage. As baselineV M depolarizes, burst frequency increases monotonically due to progressive shortening of the interburst interval and burst duration decreases monotonically due to cumulative voltage-dependent inactivation of the burst-generating currentI NaP (Fig. 2, B andC). During bursts, spike frequency is highest at burst onset then declines monotonically until burst termination (Fig.2 B), reflecting the inactivation kinetics ofI NaP.

Fig. 2.

Behaviors of the pacemaker neuron model. A: 40 s of simulated activity in the model (synaptically isolated) at severalE K, using Formula NaP = 2.8 nS and Formula L(K) = 2.4 nS. B: sample traces from A on an expanded time scale atE K = −78.9, −78, and −76.5 mV (left), and single bursts expanded and plotted side-by-side with their corresponding instantaneous spike frequency on the same time scale (middle and right).C: the frequency and duration of oscillatory bursts plotted versus E K in the cell fromA and B. E K(Min) andE K(Max) were −79 and −73.5 mV, respectively. D: burst frequency vs.E K plotted for several Formula L(K) between 1.75 and 2.75 nS ( Formula NaP = 2.4 nS). Eand F: burst frequency (E) and burst duration (F) vs. E K plotted for several Formula NaP between 2.1 and 2.9 nS ( Formula L(K) = 2.4 nS).

Maximal conductances g¯ NaPand g¯ L(K) influence the frequency and duration of bursts (Butera et al. 1999a). We define f Max andf Min as the maximum and minimum burst frequencies generated by the model,E K(Max) andE K(Min) as the highest and lowestE K that support bursting behavior, andd Max andd Min as maximum and minimum burst duration. g¯ L(K) primarily controls E K(Max) andE K(Min) but does not greatly affect burst frequency or duration. Decreasing g¯ L(K) (for fixed g¯ NaP) shifts the burst frequency-E K relationship to the right (Fig. 2 D), increasing f Max,f Min, and input resistance concomitantly. The burst duration-E Krelationship is equivalently shifted as g¯ L(K) decreases (not shown).

g¯ NaP influencesE K(Min),f Max,f Min,d Max, andd Min but notE K(Max) (Fig. 2, E andF). Increasing g¯ NaP (for fixed g¯ L(K)) supports a wider range of burst frequencies (generally higherf Max and lowerf Min) and prolongs burst duration (both d Max andd Min increase as g¯ NaP increases). Greater g¯ NaP also facilitates bursting at more hyperpolarized V M, which lowers E K(Min). However, g¯ NaP does not affectE K(Max), which is instead determined by the V M where tonic spiking emerges; spike threshold depends on I Na (notI NaP) and is therefore equivalent for all cells regardless of g¯ NaP. For all g¯ L(K) and g¯ NaP employed (Table 1) we found: f Min ≈ 0.05 Hz andf Max ≈ 0.95 Hz, andd Max ≈ 0.75 s andd Min ≈ 0.2 s (Fig. 2,E and F). The values for g¯ NaP and g¯ L(K) illustrated in Fig. 2are contained in the 95% confidence intervals for the normal distributions used to assign parameters during network simulations. Therefore cells with these representative properties participate in network simulations in this study. These properties are detailed byButera et al. (1999a).


We compared pacemaker neurons in vitro with the model. Extracellular recordings were performed to avoid altering cytosolic contents and thus intrinsic properties and to facilitate random sampling throughout the pre-BötC. Respiratory pacemakers discharged bursts of spikes coincident with XIIn motor discharge at 9 mM [K+]o (n= 28). Sixty-four percent of these neurons also discharged ectopic bursts between XIIn cycles (Fig. 3,A and B, →), as previously shown (Koshiya and Smith 1999a). Figure 3shows an inspiratory pacemaker neuron with spiking coincident with the onset of XIIn discharge (A) and a pacemaker neuron inB with preinspiratory spiking that precedes XIIn discharge (spiking precedes XIIn discharge) (Johnson et al. 1994). Spike discharge patterns were determined by inspiratory cycle-triggered spike histograms (C and D). To confirm that these cells had pacemaker properties (i.e., could burst intrinsically in the absence of rhythmic synaptic drive), we applied CNQX or low Ca2+ solution, which eliminate network activity by blocking excitatory synaptic transmission to respiratory neurons (Funk et al. 1993; Ge and Feldman 1998;Johnson et al. 1994; Koshiya and Smith 1999a). The concentration of CNQX used (20 μM) has previously been shown from whole cell recording to completely block rhythmic excitatory synaptic drive currents in pacemaker neurons (Koshiya and Smith 1999a; Del Negro, unpublished observations). CNQX or low Ca2+ conditions were maintained afterward to examine intrinsic cellular properties in the absence of respiratory rhythm.

Fig. 3.

Extracellular recordings of pacemaker neurons in vitro.A: inspiratory pacemaker neuron before (top and middle) and after (bottom) 20-μM CNQX application, which blocks rhythmic activity in the slice. B: preinspiratory pacemaker neuron before (top) and after (bottom) application of low-Ca2+ solution to block network activity.Middle: an expanded burst from above to accentuate the preinspiratory spike discharge. In both A andB, [K+]o = 9 mM, and inspiratory motor output is displayed via the integrated XIIn activity (∫XIIn). Ectopic bursts that occurred between ∫XIIn bursts are indicated (→). C and D: inspiratory cycle-triggered spike histograms (bin size =20 ms) for the inspiratory (A) and preinspiratory (B) neurons above, displayed with the synchronized cycle-triggered averages of XIIn discharge. Spike frequency histograms and XIIn averages were computed from 10 min of continuous network activity.

Neuronal excitability was manipulated in vitro to compare pacemaker cell burst frequency and duration with model predictions. Elevation of [K+]o caused quiescent neurons to begin bursting and subsequently increased burst frequency and decreased burst duration (Fig. 4). Tonic spiking occurred at the highest [K+]o (e.g., 14 mM in Fig. 4 A). This voltage-dependent behavior is consistent with the model and published data (Koshiya and Smith 1999a; Smith et al. 1991). In all cells tested, the intraburst spike frequency declined from burst onset to termination similar to the model cell (compare Fig. 4 B with 2B).

Fig. 4.

Properties of pacemaker neurons in vitro. A: 40 s of recorded activity in a representative pacemaker neuron at several [K+]o in the presence of 20-μM CNQX.B: traces from A at [K+]o = 8, 10, and 12 mM and single bursts on an expanded 600-ms time scale. Expanded bursts are displayed side-by-side with corresponding instantaneous spike frequency (○, plotted left ordinate: 0–75 Hz) and spike-frequency-time histograms (░, right ordinate, bin size = 20 ms) for all the bursts recorded during 10 min of continuous acquisition. C: burst frequency vs. [K+]o for 5 pacemaker neurons in 20 μM CNQX conditions demonstrating heterogeneity. D: mean burst duration (±SD) in sample pacemaker neurons from C.● in C and D illustrate burst frequency and duration for the cell in A and B.

In 2 of the 28 pacemaker cells sampled, we observed small spike-like discharges superimposed during bursts that were CNQX resistant. These small spikes could be due to electrotonic coupling with a neighboring respiratory neuron (Rekling et al. 2000) or voltage transients in the axon hillock or dendritic tree of the recorded neuron (Fig. 4, A and B).

In neurons for which complete data sets were obtained from quiescence to bursting and ultimately to tonic spiking (e.g., Fig. 4 A), we found these bursting characteristics:f Min ≈ 0.05 Hz andf Max ≈ 0.7 Hz, andd Max ≈ 0.6 s andd Min ≈ 0.1 s (Fig. 4,C and D, n = 7). Individual neurons exhibited heterogeneity in the range of burst frequencies and duration, i.e., in the minimum and maximum values of [K+]o at which bursting behavior was initiated and subsequently transitioned to tonic spiking. Although [K+]o cannot be directly compared with E K in simulations (Forsythe and Redman 1988), the observed monotonic increase in burst frequency and decrease in burst duration with [K+]o were consistent with the model incorporating heterogeneous g¯ NaP and g¯ L(K).

I NaP dominates the burst-generating mechanism of model 1. To evaluate the plausibility of this mechanism for the pacemaker cells, we examined whether intrinsic Ca2+ currents contributed to burst generation by switching from normal ACSF + CNQX to low-Ca2+solution containing CNQX in five experiments (Fig.5). CNQX application first caused a statistically significant decrease in burst duration (P< 0.05) due to loss of inspiratory synaptic drive currents (Fig.5 B). Subsequently switching to low-Ca2+ + CNQX solution did not significantly affect burst duration (P ≈ 0.5) nor the spike frequency-time relationship (Fig. 5 A), which is consistent with a Ca2+-independent burst-generating mechanism such as an I NaP-like current.

Fig. 5.

Effects of low-Ca2+ solution on bursting-pacemaker properties. A: an 8-s sample of pacemaker activity in 20 μM CNQX and 9 mM [K+]o before and after application of low-Ca2+ solution. Spike frequency-time histograms (bin size = 20 ms) are shown for both CNQX (top) and CNQX + low-Ca2+ conditions (bottom). B: category plot of mean burst duration (±SE) in control, CNQX, and CNQX + low-Ca2+conditions (n = 5); the statistically significant change in burst duration is shown (*, P < 0.05).

Effects of excitatory synaptic coupling on pacemaker neuron activity

In network simulations the synaptic drive currentI syn(e) generally prolongs pacemaker burst duration and decreases burst frequency compared with intrinsic cell activity in the absence ofI syn(e) (Butera et al. 1999b). To examine the role of excitatory synaptic input in vitro and evaluate these model predictions, we measured burst frequency and duration in 19 pacemaker cells before and after blocking network activity with CNQX and in 14 pacemaker neurons after blocking all chemical synaptic transmission with low-Ca2+solution.

To directly compare these experiments with model predictions, we performed new simulations using the 50-cell pacemaker network with the parameter distributions specified in Table 1. Figure6 A shows synchronized rhythmic activity in the synaptically coupled network and after removal of coupling. I syn(e) synchronizes neuronal activity to produce population-level bursts. After uncoupling, cells desynchronize their bursting or show tonic spiking or quiescence (due to heterogeneity). One pacemaker neuron exhibiting bursting in both coupled and uncoupled conditions was randomly selected in each of 19 network simulations to mimic sampling one pacemaker neuron per slice preparation in vitro. Burst frequency and duration in the sample cell were computed for coupled and uncoupled states.

Fig. 6.

Effects of removing excitatory synaptic input in model pacemaker neurons and in pacemaker neurons in vitro. A: 20 s of pacemaker-network activity during synaptically coupled and uncoupled states in the model. Top traces are running spike frequency-time histograms, computed from spike times in the 50-cell population (bin size = 10 ms). Below are raster plots for spike times in all 50 cells, sorted on the ordinate axis by cell index number. B: scatter plots of burst frequency (in Hz) before and after synaptic uncoupling, plotted on the abscissa and ordinate respectively. From left to right, the synaptic uncoupling experiments were performed in the model (n = 19), in the slice preparation using 20 μM CNQX (n = 19), and in the slice using low-Ca2+ solution (n = 14). C: scatter plots of burst duration (in s) before and after synaptic uncoupling, plotted on the abscissa and ordinate respectively. Left andmiddle: the synaptic uncoupling experiments were performed in the model (n = 19) and in the slice using 20 μM CNQX (n = 19). Right: the duration of inspiratory bursts are plotted vs. ectopic bursts in the same cell (n = 11).


Figure 6 B (middle) shows the burst frequency of pacemaker neurons in vitro before and after CNQX application on the abscissa and ordinate, respectively. From control to CNQX conditions, burst frequency decreased in the majority of cells, resulting in a significant mean decrease from 0.22 to 0.18 Hz (P < 0.05, n = 19). Contrary to these in vitro results,Butera et al. (1999b) reported that removal ofI syn(e) generally increased burst frequency in simulations. What could cause this difference between model predictions and experimental measurements? The original synaptic uncoupling simulations by Butera et al. eliminatedI syn(e) but notI tonic(e), which regulatesV M and thus voltage-dependent bursting behavior. However, pre-BötC neurons putatively receive both CNQX-sensitive tonic and phasic excitatory synaptic input (Funk et al. 1993; Ge and Feldman 1998), modeled by I tonic(e) andI syn(e), respectively. Pharmacological removal of both types of excitatory input by CNQX in vitro could explain the disparity between experiment and theory. To explore this possibility, we performed new synaptic uncoupling simulations and simultaneously blocked I syn(e) andI tonic(e) in the model to mimic the effects of CNQX (e.g., Fig. 6 A). Mean burst frequency in the 19 randomly selected pacemaker neurons decreased after blockade ofI syn(e) andI tonic(e), from 0.34 to 0.26 Hz (Fig.6 B, left, P < 0.05), which matched experimental results.

The preceding analysis assumes that the majority of synaptic inputs to pacemaker cells in vitro are excitatory, but synaptic inhibition or other endogenously active synaptic inputs that modulate pacemaker cell excitability may also influence rhythm generation. Johnson et al. (1994) used low-Ca2+ solution to block synaptic activity, and Butera et al. (1999b) used their data to show that synaptic uncoupling via low Ca2+ increased burst frequency in the majority of pacemaker neurons sampled (although statistical significance was not evaluated). These results contradict the present data using CNQX, so we postulated that there may be other sources of synaptic transmission, particularly tonic inhibition, that are blocked by low-Ca2+ solution but not by CNQX. Consistent with this idea and Johnson et al. (1994)'s original results, new data were obtained in this study using low-Ca2+ solution to synaptically uncouple pacemaker cells, which showed a mean increase in burst frequency from 0.21 to 0.28 Hz (P < 0.05, n = 14, Fig. 6 B, right). These data are consistent with the presence of CNQX-insensitive sources of synaptic transmission (acting simultaneously with tonic excitation) that ordinarily suppress burst frequency in pre-BötC pacemaker cells.


I syn(e) generally increased burst duration in model cells by contributing additional inward current during the inspiratory phase (Butera et al. 1999b). Figure 6 C (left) shows burst duration for new simulations before and after removingI tonic(e) andI syn(e) (same cells as Fig.6 B). Burst duration was generally lower after excitatory synaptic blockade: the mean decreased from 0.75 to 0.58 s (P = 0.06). This result was not statistically significant at P < 0.05 because a subset of cells showed the opposite effect, an increase in burst duration after blocking I tonic(e) andI syn(e). This subset of neurons expressed relatively large g¯ NaP (due to random assignment) and was entrained by network rhythms at burst periods too short to allow for complete time-dependent deinactivation ofI NaP during the interburst interval. We did not observe a similar subset of neurons experimentally. In vitro, burst duration consistently decreased after CNQX application; the mean decreased from 0.37 to 0.27 s (P < 0.001, n = 19; Figs. 5 B and 6 C, middle). To further examine the relationship betweenI syn(e) and burst duration in vitro, we compared inspiratory bursts to ectopic bursts that occur when the XIIn is silent. Inspiratory burst duration always exceeded ectopic burst duration, 0.37 versus 0.2 s (P < 0.001,n = 11; Fig. 6 C, right). In network simulations, highly depolarized neurons also exhibited ectopic bursts that were shorter in duration and contained fewer spikes than inspiratory bursts (Figs. 6 A, raster plots, and 7, cell 28).

Fig. 7.

Simulations of the pacemaker-network coupled to a follower population.A: 12 s of simulated activity in the pacemaker-network displayed as a running spike frequency-time histogram (bin size = 10 ms) computed from spike times in the 50-cell network (top), to show population activity (similar to ∫pre-BötC in vitro). Membrane potential trajectories are shown below population activity for cells 1, 6, 28, and45 in the pacemaker network; the voltage calibration forcell 6 applies to all pacemaker and follower cells illustrated. Below the sample voltage traces is raster plot of spike activity in all 50 cells of the pacemaker network.B: activity in the follower population (bottom) displayed as a running spike frequency-time histogram (bin size = 10 ms), which approximately mimics the ∫XIIn in vitro. Sample voltage traces for cells 1, 3, and19 are shown and a raster plot of spike activity for all 50 neurons of the follower population. Note that intrinsic heterogeneity causes dispersion of cellular spiking activity in both the pacemaker and follower populations, but the dispersion is more pronounced in pacemaker cells.

As illustrated in Fig. 6 C, the blockade of CNQX-sensitive synapses significantly decreases burst duration. After CNQX, exposing neurons to low Ca2+ + CNQX did not influence burst duration (Fig. 5 B). These data suggest that inspiratory synaptic drive currents are predominantly excitatory (CNQX-sensitive), Ca2+ currents do not significantly contribute to inspiratory burst generation, and CNQX-insensitive synapses (which are blocked by low Ca2+ and do not affect burst duration) mainly influence burst frequency.

Evaluation of population-level activity in the pre-BötC and follower neurons

We tested model predictions regarding population-level activity in the rhythm-generating kernel and in a hypothetical population of follower neurons. The pacemaker network provided excitatory synaptic drive to follower cells that attempt to model cells that transmit inspiratory drive to hypoglossal motoneurons (see methods) (see also Wilson et al. 1999). Network activity in the pacemaker and follower cells is displayed as a running histogram of action potentials in the population and as a raster plot of spike times for the two groups of 50 cells (Fig. 7). Intrinsic heterogeneity in the pacemaker population causes variability in spiking behavior: cells with low g¯ NaP and/or high g¯ L(K) show weak inspiratory activity (e.g., Fig. 7, cell 6); cells with high g¯ NaP and/or low g¯ L(K) exhibit strong inspiratory activity (cell 1) or strong inspiratory activity and ectopic bursts (cell 28). Cells with very low g¯ L(K) are highly depolarized and show inspiratory-modulated tonic spiking (cell 45).

To test model predictions regarding population-level activity, we recorded inspiratory discharge in vitro from the XIIn and locally in the pre-BötC using suction electrodes applied to the caudal surface of the slice; this exposes pacemaker neurons through the caudal boundary of the pre-BötC (Fig. 1 A). The temporal relationship of these signals was obtained from cycle-triggered averaging. Pre-BötC activity preceded XIIn activity by 100–400 ms (Fig. 8 A), consistent with the proposal that the rhythm originates in the pre-BötC (Rekling and Feldman 1998; Smith et al. 1991). Neural activity in the pre-BötC outlasted the XIIn activity for 100–200 ms (Fig. 8 A). Bilateral pre-BötC recordings were synchronous and nearly identical (Fig. 8 C).

Fig. 8.

A: comparisons of local activity patterns recorded from the pre-BötC (gray traces) and XIIn (black traces) in a typical slice preparation at [K+]o = 8, 10, and 12 mM. B: comparison of local activity patterns (running spike frequency-time histograms with bin size = 10 ms) in the model pacemaker (gray trace) and follower populations (black traces) atE K = −82 mV, analogous to superimposed plots of local ∫pre-BötC and ∫XIIn recordings in slices.C: local ∫pre-BötC recordings obtained bilaterally at 10-mM [K+]o from the slice preparation in A. Right (gray) and left (black) recordings are superimposed using normalized ordinate axes. Cycle-triggered averaging was applied to 30 consecutive bursts in the pre-BötC and XIIn in vitro to generate average traces inA and C. Inspiratory-like bursts displayed in B for the model populations were cycle-trigger averaged for 90 s of simulated time.

Analogous temporal relationships were observed in the model where we compared activity in the pacemaker population to follower cells at various E K. Rhythmic bursts in the follower population display rapid onset followed by a ramp-like activity decline similar to XIIn discharge, whereas the pacemaker population discharge both preceded and then outlasted follower activity by 100–200 ms (Fig. 8 B). The temporal dispersion of spiking activity in the pacemaker population results from heterogeneity (Butera et al. 1999b). There was less dispersion in the follower cells because they lack I NaPand thus are more homogenous (Fig. 7).


Butera et al. (1999b) found that elevating excitability in the network depresses inspiratory burst amplitude even though burst frequency concomitantly increases. This is a unique feature of the network composed of model 1 pacemaker cells and does not apply to identical networks of Butera et al. (1999a)'s model 2 neurons (data not shown). The depression of burst amplitude occurs because depolarization of model 1 pacemakers cumulatively inactivatesI NaP, causing fewer spikes per burst in constituent cells (Fig. 2) and consequently smaller bursts in the population (Butera et al. 1999b: Fig. 7 and pp. 405, 413, 414). Originally, Butera et al. used a parameter other thanE K to control excitability. Here we employed E K to mimic experimental manipulation of [K+]o. Pacemaker population activity is displayed in spike-time histograms and in plots of mean burst area versus E Kin Fig. 9, A and B. The burst-area E K data were pooled from 10 sets of simulations and normalized atE K = −88 mV. The results of these new simulations using E K matched the original study (Butera et al. 1999b, Fig. 7): net activity in the model-1-pacemaker population decreased as a function ofE K (Fig. 9 B, open circles). Figure 9 B also shows that generally 5–20% of the model pacemaker neurons were in their intrinsic bursting state at each level of E K (i.e., if synaptic coupling were eliminated these cells would continue bursting and the others would become silent or tonically active).

Fig. 9.

A: cycle-triggered averages of pacemaker-network activity and follower population activity at severalE K in a typical simulation. Traces are displayed using spike frequency-time histograms (bin size =10 ms) and a 2-s time scale. The spikes/bin calibration applies to all traces, referenced to baseline. B: normalized mean burst area (±SE) in the pacemaker population plotted vs.E K (left ordinate, ○). Follower population activity (not normalized) is also plotted vs.E K using arbitrary units (right ordinate, ●). Mean pacemaker and follower population activities were computed from 10 sets of simulations at allE K. The bottom bar graphshows the percentage of conditional pacemakers in the 50-cell network in their voltage-dependent bursting state at each level ofE K. C: cycle-triggered averages of ∫pre-BötC and ∫XIIn activity from 8 to 18 mM [K+]o in a typical slice preparation, plotted on a 2-s time scale. One-microvolt calibration applies to all traces referenced to baseline. D: normalized mean burst area (±SE) for ∫pre-BötC and ∫XIIn activity (○ and ●, respectively) in 14 slice preparations plotted vs. [K+]o.

We examined the relationship between neuronal excitability and population-level activity in the pre-BötC in vitro by varying [K+]o from 8 to 18 mM, collecting 10 min of data at each [K+]o and averaging the population bursts. Inspiratory activity was quantified using population burst area. Average traces from a representative experiment are shown (Fig. 9 C) with plots of mean burst area versus [K+]o (Fig.9 D). The burst-area [K+]o data were pooled for 14 slice experiments and normalized at 10-mM [K+]o since this concentration typically evoked maximum population discharge. Pre-BötC activity declined monotonically with increasing [K+]o, similar to the burst area-E Krelationship in the model (compare Fig. 9, B andD, ○).


To examine transmission properties of premotor circuits in vitro, we compared XIIn discharge to pre-BötC population activity at several [K+]o. Average traces are shown in Fig. 9 C with plots of burst area (Fig.9 D, n = 14 slices normalized at 10 mM [K+]o). Inspiratory XIIn activity (Fig. 9 D, ●) declined as a function of [K+]o, resembling the pre-BötC-[K+]oplot (○). This suggests that inspiratory activity is conveyed to hypoglossal motoneurons by a nearly linear transmission pathway at many levels of excitability. The model follower population failed to reproduce these results (see discussion). Instead, follower activity increased as a function of E K(Fig. 9 B, ●).

Contributions of synaptic coupling to network rhythm: experimental tests with split-slice preparations and modeling results with split-pacemaker networks

Excitatory coupling synchronizes pacemaker cell activity (Koshiya and Smith 1999a). Butera et al. (1999b) examined the role ofI syn(e) by varying the coupling conductance g¯ syn(e) and found two main effects. First, coupling strength and the network burst frequency were inversely related: weaker coupling results in higher frequency (and vice versa) (Butera et al. 1999b, Figs. 1, 2, and8). Second, coupling strength influences the cycle-to-cycle stability of network activity: weak coupling caused irregular activity patterns where net population activity fluctuated periodically (Butera et al. 1999b, their Fig. 4).

To evaluate these theoretical roles forI syn(e), we performed in vitro experiments and new simulations that reduce (but do not eliminate)I syn(e). Midline-crossing projections in the slice connect pacemaker cells in the pre-BötC bilaterally (Koshiya and Smith 1999a). We surgically ablated these connections by splitting the slice through the midline, which created two split slices that continued to oscillate independently (Fig. 10). Pre-BötC activity was bilaterally synchronous in intact slices (Figs. 8 C and10 A) but became independent and temporally uncorrelated in split slices (cross correlograms not shown).

Fig. 10.

Split-slice preparations in vitro and in the model. A: schematic drawings of the intact and split-slice preparations with ∫pre-BötC recordings obtained bilaterally. B: analogous states in the 50-cell pacemaker network before and after applying split-network conditions (see methods). Population activity is displayed separately for cells 1–25 andcells 26–50 using running spike frequency-time histograms (bin size = 100 ms).

To simulate the split-slice experiment, we synaptically uncoupled model pacemaker neurons 1–25 from 26–50, leaving all other parameters intact (includingg tonic(e), see methods). The separated subpopulations (each with n = 25) were monitored separately before and after applying the split-slice condition (Fig. 10 B) to mimic recording separately from left and right halves of the pre-BötC in vitro. Rhythmic activity in the network was synchronous when intact but became completely independent and uncorrelated in the split-network conditions, which resembled activity in the left and right split slices in vitro. In both modeling and experiments, splitting the slice/network resulted in lower population activity signal-to-noise ratios.


Effects on oscillation frequency.

To test the prediction that weaker coupling, induced by severing midline-crossing connections, increases the inspiratory frequency, we measured frequency in intact whole slices over a range of [K+]o and then repeated the experiment in split slices for comparison. Intact slices became rhythmically active at 5 mM [K+]o and the mean frequency increased monotonically until 16 mM, withf Min ≈ 0.05 Hz andf Max ≈ 0.5 Hz. At [K+]o ≥ 16 mM, mean frequency declined slightly (Fig. 11 A2, •, n = 14). Split slices became rhythmically active at 4 mM [K+]o, and the mean frequency increased monotonically as a function of [K+]o until 12 mM, withf Min ≈ 0.1 Hz andf Max ≈ 0.45 Hz. At [K+]o ≥ 12 mM, the split slice mean frequency declined slightly (Fig. 11 A, 1 and2, ○, n = 17). Splitting the slice resulted in a higher mean frequency at all [K+]o ≤ 12 mM due to a leftward shift of the frequency-[K+]o relationship (Fig.11 A2).

Fig. 11.

Properties of split slices and split-model networks. In vitro data are shown in A, 1–3, and model data in B, 1–3. A1 and B1: 30-s sample ∫pre-BötC traces in vitro and in the model at multiple [K+]o in vitro and multipleE K in the model. Model data are displayed using spike frequency-time histograms in cells 1–25 in the split network (bin size = 100 ms). A2 andB2: mean burst frequency vs. [K+]o/E K for intact (■) and split conditions (○) in vitro (A2, n = 14) and in the model (B2, n = 10). A3 andB3: coefficient of variation (mean/standard deviation) for burst area plotted vs. [K+]o/E K for intact (■) and split conditions (○) in vitro (A3) and in the model (B3).

In the intact model network, rhythmic activity occurred for −90 ≤ E K ≤ −74 mV, withf Min ≈ 0.1 Hz andf Max ≈ 0.9 Hz (Fig.11 B1). Mean frequency increased monotonically as a function of E K over this entire range (Fig.11 B2, •). In split slices, the mean frequency increased at E K > −90 mV, compared with the intact network, and no coherent split-network rhythms were observed at E K > −80 mV. The frequency increase in the split network resulted from a change in the slope of the frequency-E K relationship (Fig.11 B2, ○) rather than a leftward shift as observed experimentally (A2). Model data were pooled for 10 sets of intact and split-network simulations at allE K. Individually, some split-network simulations (and some split slices in vitro) exhibited a lower frequency after the split (e.g., Fig. 10 B), but the mean effect was an increase in frequency due to change in the slope of the frequency-E K relationship (model) or leftward shift of the frequency-[K+]o plot in vitro (Fig. 11, A2 and B2).

The model did not reproduce the frequency decline that occurred in whole and split slices in vitro at high [K+]o(Fig. 11 A2). Instead, the model displayed monotonic frequency-E K relationships for both whole and split networks (Fig. 11 B2) (Butera et al. 1999b, Figs. 7-9). This disparity can be explained by synaptic inputs active in vitro (which are not included in the model) that tend to depress network activity at the highest excitability levels (Johnson et al. 2000)

Irregular periodic activity states.

Weak I syn(e) in the model caused irregular rhythms (Butera et al. 1999b, their Fig. 4). We tested whether similar irregularities could arise in vitro due to reduced g¯ syn(e) in the split slices and in the split model network. Rhythms in the split slice in vitro and in the split model showed cycle-to-cycle variability (Figs. 10 and 11,A and B). Cycle-to-cycle stability was plotted as the coefficient of variation (CV) for burst area before (8, 10 mM [K+]o) and after the split (8–15 mM [K+]oFig. 11, A3 and B3). The model showed a >10-fold increase in the CV (more irregular) going from intact to split conditions at all E K tested (B3). A similar ∼10-fold increase in CV occurred in split slices in vitro, which was observed at all [K+]o tested (A3). The CV for burst period also increased after splitting the slice by a factor of ∼2 (not shown).


We performed experiments in vitro and new simulations to evaluate mathematical models of neonatal rat inspiratory pacemaker neurons and the excitatory pacemaker-network that is hypothesized to generate respiratory rhythm in the pre-BötC in vitro. Butera et al. (1999a,b) developed these models in the first two articles of this series; here we tested fundamental bases of the models and novel model predictions regarding cellular and population-level contributions to generation and control of respiratory rhythm. Can bursting properties of pre-BötC pacemaker neurons be replicated by the minimal mathematical model (model 1) incorporating a sole burst-generating current I NaP? How does excitatory coupling strength affect the network rhythm? Do excitatory synapses (tonic and phasic) influence pacemaker cell behaviors in the context of network activity? And how heterogeneous are pacemaker neurons and does heterogeneity affect the temporal evolution of the population-level inspiratory burst?

We found pre-BötC neurons that express silent, bursting, and tonic spiking states. Burst frequency increased and burst duration decreased as cellular excitability increased, consistent with model-1-like behavior. Intrinsic properties were heterogeneous, as reflected in cell-to-cell differences in the range of oscillation frequencies, and the levels of excitability at which cells undergo transition from silence to bursting and from bursting to tonic spiking. The burst-generating mechanism did not require intrinsic Ca2+ currents, consistent with anI NaP-like ionic mechanism. Fast excitatory synapses synchronized inspiratory bursting in the pre-BötC and generally prolonged burst duration. Inspiratory synaptic drive reduced population burst frequency, although this effect was difficult to distinguish due to multiple ongoing synaptic influences (not all of which were included in the model). Elevating neuronal excitability via [K+]o increased single cell and population burst frequency and reduced pacemaker cell and population burst duration and amplitude, respectively. Cellular heterogeneity caused temporal dispersion of activity in the rhythm-generating model network, which was mirrored in vitro by inspiratory population bursts in the pre-BötC that arise prior to (and outlast) inspiratory XIIn motor discharge. Finally, weak excitatory coupling caused irregular network rhythms.

Model and data limitations


Limitations of the mathematical models are covered in the first two articles (Butera et al. 1999a,b). Model 1 incorporates the minimum set of conductances to generate bursting similar to pacemaker cells in the pre-BötC. Real pacemaker neurons most likely express a greater number of intrinsic membrane conductances than are represented in model 1, which was developed from limited data. However, the goal was to represent essential burst-generating mechanisms and keep model 1 as simple as possible for computational tractability during network simulations. Alternative models for pacemaker cells were evaluated in the first article of the series, which favored model 1 (see Butera et al. 1999a). In this study, we more thoroughly evaluated model 1 by varying the main parameter values ( g¯ NaP and g¯ L(K)) and comparing model behavior to cells in vitro (see following text).

In constructing networks of model 1 cells, we have assumed all-to-all synaptic coupling and applied a normal distribution of coupling strengths. The actual connectivity pattern and range of coupling strengths in the neonatal rat pre-BötC is unknown, although a recent report makes quantitative predictions regarding chemical and electrotonic synapses in the neonatal mouse pre-BötC (Rekling et al. 2000). Our coupling scheme is a starting point for the neonatal rat pre-BötC that gives a low CV for respiratory network output (frequency and amplitude) similar to the intact slice at the low excitability levels. We have tested sparsity of connections by systematically varying cell-to-cell connection probability from 100 to 10% (while balancing for coupling strength), which negligibly affected network output frequency and mean inspiratory amplitude unless connection probability was <30%. Lower connection probabilities increased the CV of network inspiratory amplitude. The network currently incorporates only fast excitatory synapses; extensions of the basic model will be required to address synaptic inhibition, gap junctions, and neuromodulation (Gray et al. 1999;Johnson et al. 1996; Rekling et al. 1996b,2000; Richter et al. 1996; Smith et al. 1995, 2000). In particular, 2 of 28 pacemaker cells in this study showed some evidence consistent with gap junctional coupling after block of excitatory synaptic connections (e.g., Fig.4 B). Nevertheless their cellular behaviors were consistent with single-cell model 1 predictions. The relative importance of gap junctional coupling versus chemical synaptic coupling in the network remains to be determined.

We distributed key parameters g¯ NaP and g¯ L(K) such that constituent neurons in the network expressed a range of intrinsic behaviors in the synaptically uncoupled state (quiescence, bursting, or tonic spiking). We do not know whether our distributions for g¯ NaP and g¯ L(K) mirror the real system, but cells in vitro did exhibit similar burst frequency and duration as model cells with the range of g¯ NaP and g¯ L(K) we employed. We chose the spread of our distributions based on the following: in general too little heterogeneity reduced the frequency range of output in the coupled network, whereas too much tended to desynchronize the cellular oscillations (see also Butera et al. 1999a).


We propose that inspiratory bursting neurons coupled by fast excitatory synapses generate respiratory rhythm in vitro based on these observations: the rhythm persists after attenuation of synaptic inhibition in vitro (Feldman and Smith 1989; Gray et al. 1999; Onimaru et al. 1990); pre-BötC pacemaker neurons have been identified in neonatal rats (Johnson et al. 1994; Koshiya and Smith 1999a) and mice (Lieske et al. 2000;Thoby-Brisson and Ramirez 2000; Thoby-Brisson et al. 2000); and ionotropic excitatory synapses are required for rhythm generation (Funk et al. 1993; Ge and Feldman 1998; Smith et al. 1991) and synchronization of pacemaker cell bursting (Koshiya and Smith 1999a).

We recorded pacemaker neurons extracellularly to facilitate random sampling throughout the pre-BötC, which is unfeasible with visualized patch clamping where sampling is confined to superficial neurons (Edwards et al. 1989). Also, it allowed us to record cells for extended periods (>4 h) without rundown and loss of recording quality due to intracellular dialysis. However, this approach could not directly examine intracellular membrane potential changes or intrinsic currents, which are data ultimately needed to fully evaluate model 1.

To analyze population activity, we assumed that the measurements obtained locally with suction electrodes in the pre-BötC of slice preparations reflected the behavior of rhythmogenic neurons. We cut thin slices so that the recorded activity would measure discharge in the pre-BötC and not activity in neighboring regions of the ventral respiratory column outside of the pre-BötC. Nonpre-BötC regions of the respiratory column are functionally and anatomically distinct (Gray et al. 1999), yet still display inspiratory activity thus potentially confounding signals if one intends to measure pre-BötC neurons specifically. Even though we maximally isolated the pre-BötC for population recordings, there is a mixture of cells in this region (Johnson et al. 1994) some of which may not be rhythm generators.

In vitro preparations show slower respiratory-related rhythms compared with breathing frequency in vivo, due primarily to temperature changes and vagotomy. These issues were originally addressed by Smith et al. (1990). The kinetic parameters of the models were accordingly scaled to match the lower frequencies in vitro. We have not attempted to conduct experiments at higher temperatures and rescale the model kinetic parameters so that frequencies would be closer to those in vivo because of the lower viability of the slice preparations at elevated temperatures. This would have precluded the type of long duration experiments conducted in our studies for the detailed data-model comparisons.

Data-model comparisons: single pacemaker cell activity


Model 1 proposes a voltage-dependent bursting mechanism: transitions occur from quiescence to bursting and from bursting to tonic spiking as cells depolarize as shown from intracellular recordings (Koshiya and Smith 1999a; Smith et al. 1991;Thoby-Brisson and Ramirez 2000). Here we manipulated [K+]o in vitro to depolarize E K and thus baselineV M and test for voltage dependence. Although changes in [K+]ocannot be directly equated to E K in the model (due to non-Nernstian behavior) (Forsythe and Redman 1988), raising [K+]o in vitro andE K in the model caused similar effects: in real and model neurons bursting commenced with a frequency of ∼0.05 Hz at low [K+]o/E Kand exhibited maximum burst frequency of ∼0.9 Hz at high [K+]o/E K.

As burst frequency increased monotonically with excitability, burst duration concomitantly decreased in pacemaker neurons from ∼0.7 to ∼0.1 s, close to model predictions for the specified range of heterogeneity ( g¯ NaP and g¯ L(K)). This decrease in burst duration is not a general property of bursting models with type I topology (Bertram et al. 1995), which suggests that the proposed burst-generating mechanism for model 1 may be qualitatively correct. In support of this mechanism, slow inactivation ofI NaP in neonatal rat pacemaker neurons has been reported (Koshiya and Smith 1999b; Smith et al. 2000) and low Ca2+ solution did not prevent bursting nor diminish burst duration (Figs. 3 and 5) (see also Johnson et al. 1994). Ca2+currents are present in the pacemaker neurons (and other pre-BötC neurons) in neonatal rats (Koshiya and Smith 1999a) and mice (Elsen and Ramirez 1998) but according to our data are not essential for bursting in rat. The evidence therefore favors a burst-generating mechanism dominated byI NaP. WhetherI NaP inactivation is primarily responsible for burst termination remains to be determined, although we note that an alternative model with a different burst termination mechanism (model 2) failed to reproduce the voltage-dependent effects on burst duration (Butera et al. 1999b).

Burst duration was always slightly longer in the model compared with neurons in vitro. This is a numerical effect since in noise-free simulations the last few action potentials of a burst evolve with very long interspike intervals (Guckenheimer et al. 1997;Nayfeh and Balachandran 1995, p. 390–402). In real biological neurons, noise-based fluctuations interrupt these long-interval spikes and prematurely terminate bursts compared with corresponding model bursts. The effect is more apparent during short-duration bursts consisting of few spikes. It is also possible that we slightly overestimated the magnitude of g¯ NaP.

Model 1 behavior depends on g¯ NaP and g¯ L(K). Our parameter distributions generated activity patterns consistent with the recordings from neonatal rats. These parameters may require adjustment (or additional intrinsic currents may need to be included) for the model to generate appropriate discharge patterns for different ionic or pharmacological milieu, or in different species such as mice. For example, inspiratory pacemaker and nonpacemaker neurons of neonatal mouse express transient outward current (I A) or hyperpolarization-activated mixed-cationic current (I h) (Lieske et al. 2000; Rekling et al. 1996a; Thoby-Brisson and Ramirez 2000;Thoby-Brisson et al. 2000) which could be incorporated in model 1. More complex pacemaker cell models incorporating these currents (Smith 1996; Smith et al. 1995) can exhibit voltage-dependent behavior qualitatively similar to model 1.


Ionotropic excitatory synaptic transmission is required for rhythm generation in vitro (Funk et al. 1993; Ge and Feldman 1998). To distinguish discrete roles for excitatory synapses, Butera et al. (1999b) separately manipulated excitatory coupling strength between rhythm-generating neurons (which increased burst duration and decreased burst frequency) and tonic excitatory drive, which predominantly regulated frequency via baselineV M. These findings predicted that excitatory synaptic blockade would decrease burst duration and increase burst frequency in pacemaker neurons. This prediction required that tonic excitation, which could partially derive from pharmacologically separate receptors, remain intact.

In agreement with the first prediction, burst duration decreased in all cells tested after CNQX application presumably by eliminating depolarizing inspiratory drive currents. As stated in the preceding text, synaptic uncoupling in the model generally decreased burst duration. However, the opposite effect, an increase in burst duration after uncoupling, did occur in a subset of model cells (data above the diagonal in Fig. 6 B). This seemingly contradictory effect arose from network entrainment at a frequency that prevented full de-inactivation of I NaP during the interburst interval. After uncoupling from the network, these cells burst at lower intrinsic frequencies, which extended the recovery time for I NaP between burst cycles and thus augmented the available I NaP and prolonged burst duration. We did not observe a subset of neurons that increased burst duration after CNQX in vitro. This disparity could have resulted from sampling bias since inspiratory neurons with strong bursts at 9-mM [K+]o were generally selected during experiments; model neurons with strong bursts in control conditions invariably decreased their burst duration after synaptic uncoupling (e.g., data points in the bottom right portion of Fig. 6 B).

Both model and real pacemaker neurons generated ectopic bursts during the interval between synchronized inspiratory bursts. In the model, synaptic input during the ectopic burst is weak or absent compared with inspiratory synaptic drive, thus the majority of the ectopic burst-generating current is intrinsic (I NaP). Inspiratory burst duration in real pacemakers always exceeded ectopic burst duration, which we attribute (based on the model) to strong synaptic drive during inspiratory bursts and a predominantly intrinsic mechanism for ectopic burst generation. All of the ectopic burst-generating cells in vitro were confirmed to be pacemakers after applying CNQX or low-Ca2+ solution. Therefore ectopic bursts can be used to identify pacemaker neurons in the functionally intact respiratory network in vitro.

The effect of removing inspiratory synaptic input on burst frequency varied: burst frequency increased in some cells and decreased in others. The mean effect was a decrease in frequency, the opposite of what was predicted. The model generated similar results when we eliminated both I tonic(e) andI syn(e). This suggests that CNQX application in vitro eliminated both phasic (inspiratory) and tonic sources of excitatory input, which is consistent with evidence that non-NMDA receptors are recruited during the inspiratory phase and may also be tonically active during rhythm generation (Ge and Feldman 1998). We conclude that inspiratory frequency in vitro may not be controlled primarily by synaptic coupling strength as predicted by the network model at medium-high coupling strengths. Nevertheless these data are consistent with the network model when coupling strength is relatively low (Butera et al. 1999b). Nevertheless like in isolated pacemaker neurons, the network frequency does depend on cellular factors that affect baselineV M such as tonic excitation,E K, or neuromessengers (e.g.,Gray et al. 1999; Johnson et al. 1996;Rekling et al. 1996b; Richter et al. 1997).

When low-Ca2+ solution was used to uncouple pacemaker neurons burst frequency generally increased, unlike the effect of CNQX in vitro or removal of bothI tonic(e) andI syn(e) in the model. Because low-Ca2+ solution blocks all forms of chemical synaptic transmission, these experiments suggest that CNQX-insensitive synaptic receptors are also activated endogenously in vitro, and the net result is to normally depress burst frequency. Some inhibitory ionotropic receptors that are endogenously activated have been identified (Johnson et al. 2001).

Population-level activity in the pre-BötC


We showed that pre-BötC activity develops prior to the onset of respiratory-related motor activity recorded from the XIIn, which is required if the respiratory rhythm originates in the pre-BötC (Rekling and Feldman 1998;Smith et al. 1991). The envelope of pre-BötC activity was bell-shaped in contrast to the XIIn pattern, which displays a more rapid onset followed by a ramp-like decline to baseline (originally analyzed by Smith et al. 1990). Activity in the pre-BötC also outlasted XIIn discharge for 100–200 ms. These temporal relationships were observed at all [K+]otested in vitro and were predicted by the pacemaker network coupled to a follower population.

The envelope of population activity is determined in the model by pacemaker neurons that possess I NaPand are heterogeneous, which causes temporal dispersion in population activity. In contrast, follower neurons are hypothesized to partially comprise rhythmic drive-transmission circuits that cannot generate bursts intrinsically and have less intrinsic heterogeneity because they lack I NaP. Therefore followers burst only in response to synaptic input resulting in a more abrupt onset of action potential discharge during the inspiratory phase (determined by synaptic not intrinsic currents) and an earlier decline in the inspiratory discharge compared with pacemaker neurons with sufficientI NaP to boost inspiratory discharge. Our follower neuron population model is hypothetical, but given the similarity of follower population discharge and activity patterns recorded from the XIIn in vitro, we postulate that neurons of the rhythmic drive-transmission circuits probably do not express intrinsic burst-generating currents.


Butera et al. (1999b) predicted that burst amplitude in a network of model 1 neurons declines monotonically as neuronal excitability increases; this does not apply to a network of model 2 pacemakers nor other generic pacemaker models (unpublished simulations). Here we performed new simulations usingE K to compare network simulations with experiments and obtained equivalent results. Our tests affirmed the model predictions: at [K+]o, >10 mM the amplitude of pre-BötC population activity declined monotonically as a function of [K+]o. In simulations, this decline in population activity relates directly to cumulative voltage-dependent inactivation ofI NaP, the same voltage-dependent mechanism that decreases burst duration in isolated model 1 pacemaker neurons. Even though no more than 30% (usually 5–20%, Fig.9 B) of the model neurons are in their voltage-dependent bursting state at a given level of E K, the net population activity nevertheless reflects the voltage-dependent inactivation properties of I NaP. This suggests that the burst-generating mechanism may be qualitatively and/or partially correct.


A fundamental prediction of the model is that tonic excitation of the pacemaker cells will increase the network oscillation frequency monotonically, a principle we tested in vitro by manipulating [K+]o to depolarize pacemaker neurons. In general, frequency increased with elevated [K+]o except at high levels of [K+]o/E K. At 17–18 mM in intact slices (and at 14–15 mM [K+]o in split slices), the inspiratory frequency declined slightly, resulting in nonmonotonic frequency versus [K+]orelationships. More recent experimental results show the predicted monotonic relationship when population activity is recorded in pre-BötC “islands” isolated from the thin slices (Johnson et al. 2001). We therefore attribute the depression of respiratory frequency at high [K+]o to synaptic inputs arising from cells extrinsic to the pre-BötC in the slice that tend to depress excitability at high [K+]o. Bicuculline-sensitive GABAA receptors appear to be activated in the pre-BötC and affect the level of excitability at which the pacemaker cells burst (Johnson et al. 2001). Thus it is clear that nonexcitatory synaptic inputs to pacemaker cells, which have not yet been included in the model, are endogenously active in vitro and cause the deviations from model predictions.


The slice was bisected experimentally to test specific predictions that a reduction, but not elimination, ofI syn(e) would increase respiratory frequency and decrease the regularity of rhythmic output (Butera et al. 1999b). Rhythm was bilaterally synchronous in the pre-BötC until splitting the slice, when it became completely uncorrelated presumably via ablation of synchronizing connections. Therefore severing the midline likely reduced coupling strength as we intended and reduced the rhythm-generating population by approximately half. Based on the additional effect of changing population size, we simulated the split-network model to examine whether the original predictions regarding reductions in coupling strength (Butera et al. 1999b) applied if population size and coupling strength decreased simultaneously. Split-slice simulations confirmed both effects: an increase in frequency and a decrease in the regularity of the rhythm.

Split-slice experiments in vitro showed similar, but not identical, effects. After midline severance, the rhythm recorded locally in the pre-BötC normally developed cycle-to-cycle variability. Weakened synaptic coupling is insufficient to completely synchronize model pacemaker cells, which possess autonomous rhythmicity and are heterogeneous. Even with normal coupling strength in the model, subsets of pacemaker neurons developed partially synchronized ectopic bursts (e.g., Fig. 7), which shows their predilection for intrinsic rhythms. When coupling is relatively weak, like in the split slice, subsets of cells can oscillate semi-independently from the fundamental network rhythm, as shown in the original paper (Butera et al. 1999b, their Fig. 4).

Cycle-to-cycle fluctuations in the population-level burst can be attributed to reductions in coupling strength. However, reducing the size of the population (as also occurs in a split slice) also promotes desynchronization of pacemaker neuron activity. Any cell whose intrinsic behavior diverges from the rest of the network becomes more apparent (and more influential) in a population half the original size.

Inspiratory frequency also increased as a result of splitting the slice in vitro. However, this effect probably does not result directly from the predicted mechanism (i.e., a decrease in excitatory coupling strength among rhythm-generating neurons). Splitting the network model increased the slope of the burst frequency-E K relationship, which can be attributed to coupling strength reduction since tonic excitatory drive was unchanged in the simulations. In vitro, at all [K+]o <14 mM, mean respiratory frequency increased after splitting the slice. This effect was due to a small leftward shift of the burst frequency-[K+]orelationship not a change in slope, suggesting some decrease in tonic excitation due to the split that caused a shift in frequency at all [K+]o rather than a main frequency effect due to reducedI syn(e).


Follower neurons were included in the model to begin constructing rhythmic drive-transmission circuits, which convey inspiratory rhythm from pre-BötC neurons to respiratory hypoglossal motoneurons. The membrane properties of these cells are unknown, although excitatory amino-acid receptors are involved in inspiratory drive transmission (Funk et al. 1993; Ge and Feldman 1998). Therefore we parsimoniously constructed the follower population using model 1 neurons with g¯ NaP = 0 nS (which precludes bursting activity) and with fast excitatory synapses as employed in the pacemaker model. If follower neurons hadI NaP, they would be identical to “weak” pacemaker cells in the kernel and thus not a distinct population.

A major disparity between the model follower population and the XIIn activity in vitro emerged as we variedE K. XIIn discharge declined monotonically with increasing [K+]o, similar to pre-BötC activity. In contrast, follower population activity increased as E K depolarized the model follower cells. This follower behavior resulted from the lack of I NaP. Intrinsic currents with voltage-dependent inactivation properties (such asI NaP) tend to depress neuronal discharge as E K increases. We posit that interneurons of the rhythmic drive-transmission circuits in vitro express some intrinsic regulatory properties that decrease their discharge activity as [K+]o increases during experiments. This mechanism could be intrinsic or synaptic in nature but probably is not a burst-generating current with voltage-dependence and kinetics like I NaP otherwise the temporal pattern of activity recorded from the XIIn would resemble the bell-shaped pattern of pre-BötC.


Our tests of model predictions affirm that the burst-generating mechanism in inspiratory pacemaker neurons is a Ca2+-independent persistent Na+ current, that pacemaker neurons are heterogeneous, and that a heterogeneous network of pacemakers coupled by excitatory synapses can account from many aspects of the inspiratory rhythm and pattern generated in vitro. These results support many ofButera et al. (1999a,b)'s conclusions regarding pacemaker cell behavior and the roles of excitatory synapses and of cellular heterogeneity in rhythm generation in vitro.


C. A. Del Negro was supported by a National Institute of Neurological Disorders and Stroke Intramural Competitive Fellowship Award. S. M. Johnson was supported in part by Office of Naval Research Grant N00014-94-0523.


  • * S. M. Johnson contributed experimental data and R. J. Butera contributed simulations and data analyses.

  • Address for reprint requests: J. C. Smith, 49 Convent Dr., Rm. 3A50, Bethesda, MD 20892-4455 (E-mail:jsmith{at}


View Abstract