Journal of Neurophysiology

Error message

Notice: PHP Error: Undefined index: custom_texts in highwire_highwire_corrections_content_type_render() (line 33 of /opt/sites/jnl-jn/drupal-highwire/releases/20151124215058/modules/highwire/plugins/content_types/

Models of Respiratory Rhythm Generation in the Pre-Bötzinger Complex. I. Bursting Pacemaker Neurons

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


A network of oscillatory bursting neurons with excitatory coupling is hypothesized to define the primary kernel for respiratory rhythm generation in the pre-Bötzinger complex (pre-BötC) in mammals. Two minimal models of these neurons are proposed. In model 1, bursting arises via fast activation and slow inactivation of a persistent Na+ current I NaP-h. In model 2, bursting arises via a fast-activating persistent Na+ current INaP and slow activation of a K+ current IKS. In both models, action potentials are generated via fast Na+ and K+currents. The two models have few differences in parameters to facilitate a rigorous comparison of the two different burst-generating mechanisms. Both models are consistent with many of the dynamic features of electrophysiological recordings from pre-BötC oscillatory bursting neurons in vitro, including voltage-dependent activity modes (silence, bursting, and beating), a voltage-dependent burst frequency that can vary from 0.05 to >1 Hz, and a decaying spike frequency during bursting. These results are robust and persist across a wide range of parameter values for both models. However, the dynamics of model 1 are more consistent with experimental data in that the burst duration decreases as the baseline membrane potential is depolarized and the model has a relatively flat membrane potential trajectory during the interburst interval. We propose several experimental tests to demonstrate the validity of either model and to differentiate between the two mechanisms.


Breathing movements in mammals are generated by networks of neurons in the lower brain stem that produce a rhythmic pattern of neural activity. Recently the primary neuronal kernel for rhythm generation has been located in the pre-Bötzinger complex (pre-BötC), a subregion of the ventrolateral medulla (Smith et al. 1991). This discovery lead to the development of rhythmic in vitro medullary slice preparations from neonatal and juvenile rodents (Funk et al. 1994;Ramirez et al. 1996; Smith et al. 1991) that capture this kernel and have become important experimental preparations for analysis of cellular and network mechanisms of rhythm generation. Current evidence (reviewed in Rekling and Feldman 1998;Smith 1997; Smith et al. 1995) indicates that rhythm generation in these slice preparations, as well as in more en bloc in vitro preparations, arises from a population of pre-BötC excitatory interneurons with intrinsic oscillatory bursting or pacemaker-like properties. It thus has become clear that to understand respiratory rhythm generation, at least in vitro, mechanisms incorporating intrinsic cellular pacemaker properties must be analyzed. Accordingly, a new mechanistic model, the hybrid pacemaker-network model (Smith 1997; Smith et al. 1995), has been proposed in which rhythm arises from the dynamic interactions of both intrinsic and synaptic properties within a bilaterally distributed population of coupled bursting pacemaker neurons. In this and the following paper (Butera et al. 1999), we present computational versions of this “hybrid” model that provide an initial analytic framework for analyzing the potential roles of cellular and synaptic processes in the generation and control of rhythm.

The hybrid pacemaker-network model departs from previous network-based models (Balis et al. 1994; Botros and Bruce 1990; Duffin 1991; Gottschalk et al. 1994; Ogilvie et al. 1992; Rybak et al. 1997) in which respiratory rhythm has been postulated to arise mainly from network interactions, particularly inhibitory connections; synaptic interactions are proposed to operate cooperatively with intrinsic cellular properties of specific classes of interneurons to produce the phase transitions required for a network-based respiratory rhythm (see Ramirez and Richter 1996; Richter et al. 1992). In these models, rhythmicity ceases when synaptic inhibition is blocked. However, in the in vitro slice and en bloc preparations, inspiratory phase respiratory activity persists when inhibitory synaptic connections are blocked pharmacologically (Feldman and Smith 1989; Ramirez et al. 1996; Shao and Feldman 1997), and neurons with intrinsic bursting oscillations have been identified (Johnson et al. 1994; Smith et al. 1991). In the hybrid model, inhibitory interactions are not essential; thus the fundamental difference between this model and previous models is that a population of synaptically coupled excitatory pacemaker-like neurons generates the inspiratory phase of respiratory network activity. This rhythm and inspiratory burst-generating kernel is embedded in a complex network, however, that provides excitatory and inhibitory synaptic mechanisms for control of oscillatory bursting as well as for generation of the complete pattern of respiratory network activity including the phasic firing of neurons during expiration (see discussion in Smith 1997; Smith et al. 1995). Our focus in this and the following paper is on modeling the rhythm and inspiratory burst-generating kernel operating in vitro. These models serve as the basis for the development of a more complete hybrid model of the respiratory network incorporating both rhythm-generating kernel and pattern-formation networks (Smith 1997) that must be included to account for inspiratory and expiratory patterns of activity in vitro and in vivo.

There currently is limited experimental information on the ionic and synaptic mechanisms generating and controlling the rhythmic bursting of pre-BötC pacemaker cells. We therefore have formulated minimal models for the pacemaker neurons with reduced parameter sets that nevertheless retain the essence of what has been measured or hypothesized for the cellular properties and synaptic interactions. Our objective has been to formulate the models in a way that facilitates exploration of general principles and mechanisms for oscillatory burst generation. In this paper, we present our pacemaker neuron models and address questions about membrane conductance mechanisms that may underlie the voltage-dependent oscillatory behavior of the candidate pre-BötC pacemaker neurons found in vitro. We conducted a systematic analysis of potential mechanisms regulating oscillatory bursting, burst frequency, and burst duration at the single neuron level. We also have derived tests that allow several general mechanisms to be distinguished that will guide experimental measurements. In the succeeding paper (Butera et al. 1999), we extend the analysis to the cell population level and consider a population of synaptically coupled bursting pacemaker neurons—a model of the inspiratory rhythm-generating kernel in in vitro preparations. We address a number of issues about the dynamics of this pacemaker network kernel, including how bursting of neurons within the population is synchronized and how synaptic interactions and intrinsic cellular properties dynamically regulate inspiratory burst frequency and duration. Preliminary reports of these modeling results have been presented in condensed form (Butera et al. 1997b,1998a,b).


All simulations were performed on IBM RS/6000 or Pentium-based UNIX/LINUX workstations. Most simulations were coded in the C programming language using the numerical integration package CVODE (Cohen and Hindmarsh 1996) available at For final simulations, relative and absolute error tolerances were ≤10−6 for all state variables. Some simulations also were performed using the interactive differential equation simulation package XPP available at

Model development

The primary features of oscillatory bursting of the candidate inspiratory pacemaker neurons in the pre-BötC (recorded from in vitro transverse medullary slice preparations from neonatal rat) are illustrated in Fig. 1 (see alsoKoshiya and Smith 1999; Smith et al. 1991). After block of synaptic transmission by low-Ca2+ conditions in the slice bathing medium, the neuron exhibits voltage-dependent oscillatory bursting behavior. As the cell is depolarized by a steady applied current under whole cell current clamp, the cell undergoes a transition from a state of rest to a state of oscillatory bursting. As the cell is depolarized further, the burst period, as well as the burst duration (Fig. 1 B), decreases. Additional depolarization maintains the cell above the action potential threshold and causes a transition to a state of beating (i.e., tonic spiking). Another salient feature is a steadily decreasing spike frequency throughout the duration of the burst (see Fig. 1 B) (see also the spike-frequency histograms of Johnson et al. 1994). Thus these neurons are presumed to have multiple functional states (quiescence, oscillatory bursting, and beating); they have been referred to as conditional pacemaker neurons (Smith 1997; Smith et al. 1991, 1995) to indicate that particular conditions must exist (ranges of depolarization level or magnitudes of burst-generating conductances) for oscillatory pacemaker-like bursting to occur.

Fig. 1.

Example of voltage-dependent properties of pre-Bötzinger complex (pre-BötC) inspiratory bursting neurons. Traces show whole cell patch-clamp recordings from a single candidate pacemaker neuron in the pre-BötC of a 400-μm-thick neonatal rat transverse medullary slice with rhythmically active respiratory network. Recordings inA and B were obtained respectively before and after block of synaptic transmission by low Ca2+conditions identical to those described in Johnson et al. (1994) (i.e., 0.2 mM Ca2+, 4 mM Mg2+, 9 mM K+ in slice bathing solution). Patch pipette solution and procedure for whole cell recording were as described previously (Smith et al. 1991, 1992). Before block of synaptic transmission, the neuron bursts in synchrony with the inspiratory phase of network activity as monitored by the inspiratory discharge recorded on the hypoglossal (XII) nerve (Smith et al. 1991). After block of synaptic activity (30 min under low-Ca2+conditions), the cell exhibits intrinsic voltage-dependent oscillatory behavior. As the cell is depolarized by constant applied current, it undergoes a transition from silence (baseline potential below −65 mV,left) to oscillatory bursting to beating (baseline potential above −45 mV, right). In the bursting regime, the burst period and duration decreases (see expanded time-base traces in B) as the baseline membrane potential is depolarized.

In investigating the biophysical basis of these features, a fast depolarizing mechanism for burst initiation and a slower opposing mechanism for burst termination must be considered. In the two models developed in the following text, initiation occurs by the activation of a persistent Na+ current (I NaP). However, burst termination in the two models occurs by contrasting mechanisms. In model 1, burst termination occurs by the slow-inactivation of I NaP-h, a persistent Na+-current with slow inactivation. In model 2,a slowly activating K+-current (I KS) is responsible for burst termination. Results for both models will be presented in this paper. These models cover the two major mechanisms (slow inactivation of inward current, slow activation of outward current) for “type I bursting” (Bertram et al. 1995a; Rinzel and Lee 1987) to occur in oscillatory cells. Detailed justification for our choices of conductance mechanisms is presented in the discussion.


Our minimal model is based on a single-compartment Hodgkin-Huxley (HH) formalism. The model’s dynamics are described completely by an autonomous set of differential equations. The time course of the membrane potential is obtained by applying Kirchoff’s current law to a single compartment neuron. In this case, the transmembrane current is equal to the sum of the intrinsic and externally applied currents, as follows:CdVdt=INaPINaIKILItonic­e+Iapp Equation 1where C is the whole cell capacitance (pF),V is membrane potential (mV), and t is time (ms). The ionic currents on the right-hand side are described in the following text. C is set to 21 pF, consistent with whole cell capacitance measurements from inspiratory neurons in vitro (Smith et al. 1992).

The conductances of the ionic currents are regulated by voltage-dependent activation and inactivation variables. The dynamics of a gating variable x are described according todxdt=x(V)xτx(V) Equation 2 x(V)={1+exp[(Vθx)/ςx]}1 Equation 3 τx=τ̄x/cosh[(Vθx)/(2ςx)] Equation 4where x (V) is the steady-state voltage-dependent (in)activation function of xand τx(V) is the voltage-dependent time constant. x (V) is a sigmoid with a half-(in)activation at V = θx and a slope that is proportional to 1/ςx(V) is a bell-shaped curve that has a maximal value τ̄x atV = θx and a half-width determined by ςx . Thus each gating variable is described by only three parameters. The preceding formulation also may be expressed in the HH formalism dx/dt = αx(V)(1 − x) − βx(V)x, where αx(V) = [1/(2τ̄x)]exp[−(V − θx)/2ςx)] and βx(V) = [1/(2τ̄x)]exp[(V − θx)/(2ςx)]. This is the simplest possible formulation for the voltage dependency of α(V) and β(V), and corresponds to assuming a single rate-limiting energy barrier between the open and close states of the channel (Jack et al. 1975).

Action potentials in the model are generated by a fast Na+current (I Na) 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 descriptionINa=ḡNam3(V)(1n)(VENa) Equation 5 IK=ḡKn4(VEK) Equation 6where the parameters of I Na are Na = 28 nS, E Na = 50 mV, θm = −34 mV, and ςm = −5 mV and the parameters of IK are K = 11.2 nS, E K = −85 mV, θn = −29 mV, ςn = −4 mV, τ̄n = 10 ms. The formulation for I Na includes two simplifications compared with the standard HH formulation. First, it is assumed that m activates so fast that it can be considered to be instantaneous without significantly altering the dynamics of the model (Krinskii and Kokoz 1973; Rinzel 1985). Second, the time course of h is assumed to be of similar dynamics asn and is approximated by h = 1 −n (Krinskii and Kokoz 1973; Rinzel 1985). The voltage-dependent functions for m 3(V) andh (V) forI Na andn 4(V) and τn(V) for I K are illustrated in Fig. 2 A.Quantitative measurements of the gating characteristics ofI Na and I K have not been made in respiratory neurons, and we are omitting other ionic currents that may modulate spike activity. Therefore we have chosen a minimal formulation that generates spike frequencies that are consistent with experimental recordings. Our formulation forI NaP-h in model 1 is identical to that for I NaP in model 2 except thatI NaP-h also has the inactivation termh. The following discussion regarding the activation parameters of I NaP-h applies toI NaP (in model 2) as well.I NaP-h is of the formINaP­h=ḡNaPm(V)h(VENa) Equation 7where NaP = 2.8 nS, θm = −40 mV, ςm = −6 mV, θh= −48 mV, ςh = 6 mV, τ̄h =10,000 ms, and E Na = 50 mV. The activation m is approximated as instantaneous because it occurs on a time scale of several milliseconds, which is several orders of magnitude faster than the time scale of h. Our characterization of the activation (m) and inactivation (h) ofI NaP-h is distinct from that ofI NaP, where similar variable names were used. Throughout the rest of this text, any reference to hwill always refer to the inactivation of I NaP-hand not I Na.

Fig. 2.

Gating and I-V characteristics of components ofmodels 1 and 2. A: spike generating kinetics. m 3(V) andh (V) ofI Na andn (V) and τn(V) of I K: note that h = 1 − n.B1: gating characteristics of INaP:m (V),h (V), and τh(V) (bold). Left:y axis scale for steady-state gating functions.Right: y axis scale for τh(V). B2:I-V plots of I NaP forh =h (V) andh = 1. First case results in a small window current at subthreshold potentials. Second case corresponds toI NaP-h with complete removal of inactivation. C1: gating characteristics ofI KS:k (V) and τk(V) (bold). Left:y axis scale for activation function;right: y axis scale for τk(V). C2:I-V plots of I NaP +I KS for k =k (V) andk = 0. First case results in a small current at subthreshold potentials. Second case corresponds toI NaP with complete removal of the opposingI KS.

Figure 2 B1 illustrates the functionsm (V),h (V), and τh(V) for I NaP-h. Figure 2 B2 illustrates I NaP-h whenh = h (V) andh = 1. This figure predicts theI NaP-h that we expect would be elicited under voltage clamp by starting at a hyperpolarized membrane potential and applying a slow voltage ramp and a fast voltage ramp, respectively, when all other intrinsic membrane currents are eliminated. Model 1 also has a passive leakage current I L, defined as I L = L(VE L). This current is K+ dominated. The nominal parameters for I L are L = 2.8 nS andE L=−65 mV.

The pacemaker neurons in the pre-BötC are proposed to receive both excitatory and inhibitory input from populations of tonically firing cells (see discussions in Smith 1997;Smith et al. 1995). For simplicity we only consider a mean level of excitatory input, represented by the conductanceg tonic-e, although analogous results are obtained if a mean level of inhibitory input is considered as well. Synaptic input from this beating population is assumed to be mediated by non-N-methyl-d-aspartate (non-NMDA) glutamatergic receptors (Funk et al. 1993; Smith et al. 1995), modeled as I tonic-e =g tonic-e(VE syn-e), where E syn-e is 0 mV, the reversal potential of non-NMDA glutamatergic synaptic currents. Throughout this paper I tonic-e is considered to be inactive (g tonic-e = 0 nS) unless otherwise indicated. I tonic-e is intended only to represent external tonic drive and not the more complex phasic synaptic interactions of these bursting neurons within the pacemaker population (modeled in the companion paper) or other types of modulatory synaptic inputs.

The applied stimulus current I app is included for completeness and is set to 0 pA unless otherwise noted.


Because there are not sufficient quantitative data to fully constrain our parameter choices, we formulated model 2 as similar as possible to model 1, enabling a fair and complete assessment of the dynamic behavior of both models. In model 2, we introduce an additional slow K+ current,I KS. Thus the current balance equation isCdVdt=INaPIKSINaIKILItonic­e+Iapp Equation 8where C is the same as in model 1. The ionic currents I Na, I K,I L, I tonic-e, andI app are identical to those in model 1. I NaP is identical toI NaP-h in model 1 except that it does not possess the inactivation term h. IKS=ḡKSk(VEK) Equation 9where KS = 5.6 nS, θk = −38 mV, ςk = −6 mV,τ̄k = 10,000 ms, andE K = −85 mV. The voltage dependency of the gating kinetics for k was chosen so thatI KS was activated by the firing of action potentials and was also dynamically active in the subthreshold potential range from −60 to −45 mV. Figure 2 C1 illustratesk (V) and τk(V). Figure 2 C2 illustratesI NaP + I KS whenk = k (V) (slow voltage ramp with other membrane currents eliminated) andk = 0 (removal of I KS, similar to a fast voltage ramp). The shapes of these plots ofI NaP + I KS are similar to those of I NaP-h illustrated in Fig.2 B2.

The slow variables h in model 1 and kin model 2 have been modeled as voltage-dependent processes based on the existence of oscillatory bursting inspiratory neurons in a low Ca2+ medium (Johnson et al. 1994). However, it is not essential for burst generation that gating processes depend directly on membrane potential. For example, similar dynamics may be achieved from models 1 and 2 if the slow gating variables were driven by the accumulation and removal of an intracellular ion such as Ca2+, although such a mechanism is not suggested by current experimental data.


Voltage-dependent behavior

There are several methods through which the level of depolarization, and thus the model of activity exhibited by the bursting neurons, may be controlled. Some of these methods are illustrated in Fig. 3. Note that changingI app to control the neuron’s level of depolarization is equivalent to varying E L(I app = LΔE L).

Fig. 3.

Mechanisms for controlling voltage-dependent activity in a single bursting neuron. Baseline membrane potential of the model may be controlled by application of an extrinsic current pulse (I app), changing [K+]o (which in turn affectsE K) or adjusting the level of tonic synaptic input, the mean level of which is represented byg tonic-e.

Figure 4 illustrates the effect of varying E L on the oscillatory activity ofmodel 1. The time course of membrane potential is illustrated in Fig. 4 A, 1–4, and the time course of the I NaP-h inactivation (h) and instantaneous membrane conductance are illustrated in Fig. 4 B, 1–4. As E L is increased, the model makes a transition from silence (A1) to bursting (A, 2 and 3). One burst cycle consists of an active phase, denoted by the repetitive firing of action potentials, and a silent phase, where the membrane potential varies slowly at a hyperpolarized value. The duration of the active phase is referred to as the burst duration, and the duration of the silent phase is referred to as the interburst interval. The burst period is the sum of the burst duration and interburst interval. The burst period and burst duration decrease and the cell becomes less hyperpolarized during the silent phase asE L is increased. V min, the minimum membrane potential encountered during the silent phase of the burst cycle, varies from −58 mV at the onset of bursting activity (E L = −60.5 mV) to −48 mV (E L = −57 mV) at the transition between bursting and beating activity. At all levels of depolarization the spike frequency decreases throughout the burst (e.g., Fig. 4, C, 2 and 3, D, 2 and 3). As the cell model becomes further depolarized, it undergoes a transition to tonic beating behavior (A4). All of these features are similar to those shown in the experimental recording of Fig. 1. During bursting,I NaP-h is inactivated progressively with the firing of each action potential (Fig. 4D, 2 and3). During the interburst phase, inactivation is removed and the instantaneous membrane conductance gradually increases asI NaP-h recovers from inactivation. Similar results to those in Fig. 4 A are obtained ifE L is separated into K+ and Na+ components and E K is varied.

Fig. 4.

Dynamic response of model 1 as a function ofE L. Left: (A, 1–4) membrane potential. Middle: (B, 1–4) instantaneous membrane conductance (g m) and inactivation ofI NaP-h (h) as cell is depolarized (E L = −65, −60, −57.5, and −54 mV for 1–4, respectively). C2/D2and C3/D3: expanded time course of V andh and the spike-frequency profile (○) during 1 burst in A2/B2 and A3/B3.Instantaneous conductance plots are truncated during the firing of action potentials.

Modulating the level of tonic drive can bias the neurons to different activity modes. This is analyzed in our model by varyingg tonic-e, which represents the mean level of excitatory input to the bursting neuron population from a tonically beating cell population (Fig. 3). The depolarizing effects ofg tonic-e on the mode of activity and pattern of bursting are similar to those of E L (quantified in the following text).

Figure 5 illustrates the effect of varying E L on the oscillatory activity ofmodel 2. The time course of membrane potential is illustrated in Fig. 5 A, 1–4, and the time course of the I KS activation and instantaneous membrane conductance are illustrated in B, 1–4. The dynamics of model 2 are similar in many respects to model 1. As E L is increased, the model progresses through regimes of silence, bursting, and beating activity. During oscillatory bursting, the burst is initiated byI NaP and terminates whenI KS has been activated sufficiently. The spike frequency decays throughout the active phase of the burst. The currentI KS is activated slowly and progressively with the firing of each action potential and deactivates during the interburst interval. The instantaneous membrane conductance (g m) decreases throughout most of the silent phase and increases toward the end of the phase. This differs frommodel 1, where g m steadily increases throughout the silent phase. This difference in the time course ofg m may provide a means for differentiating between the mechanisms of the two models in vitro (seediscussion).

Fig. 5.

Dynamic response of model 2 as a function ofE L. Left: (A, 1–4) membrane potential. Right: (B, 1–4) instantaneous membrane conductance (g) and activation of I KS (k) as cell is depolarized (E L = −65, −59.5, −50, and −40 mV for 1–4, respectively). C2/D2and C3/D3: expanded time course of V andk and the spike-frequency profile (○) during 1 burst in A2/B2 and A3/B3.Instantaneous conductance plots are truncated during the firing of action potentials.

The dynamics of model 2 were investigated across a wide range of parameters for I KS andI NaP. For all parameter ranges where bursting activity was supported, we found three key differences from the dynamics of model 1. The membrane potential trajectory was not as flat during the interburst interval, the burst duration did not decrease with depolarization (Fig. 8), and g mdid not exhibit a monotonic trend during the silent phase.

Robustness of frequency control

Models 1 and 2 support burst frequencies that vary over at least an order of magnitude withE L or g tonic-e, and this range of burst frequencies is robust across a range of parameter values. Figure 6 quantifies these effects for model 1 for various values of NaP. The burst duration, burst period, and frequency of bursting and beating are plotted asE L (Fig. 6 A) andg tonic-e (Fig. 6 B) are varied. Equivalent values of I app withE L fixed to −65 mV are also shown onA.

Fig. 6.

Intrinsic and extrinsic mechanisms for frequency control inmodel 1. Panels show effect of varyingE L (A) andg tonic-e (B) on burst duration, burst period, burst frequency, and beat frequency. These curves were calculated for multiple values of NaP. Bursting activity does not occur for NaP = 2.0 nS in A and for NaP = 2.0 and 2.4 nS in B. I APP below E L scale shows equivalent variation in I APP withE L fixed to −65 mV.

In both panels, there is a minimum value of NaP for which oscillatory bursting can occur. If NaP is too low, the only supported modes of activity are quiescence and beating as shown formodel 1 in Fig. 7. For all values of g NaP where oscillatory bursting occurs, the burst frequency may be varied by approximately an order of magnitude, and the burst period and burst duration decrease as the neuron is depolarized. This is consistent with the data illustrated in Fig. 1. Increasing NaP increases the dynamic range of values of E L where bursting is supported. We also have observed (not shown) that increasing NaP increases the range of values ofV min encountered as a function ofE L during oscillatory bursting. At a given level of E L, increasing NaP increases burst frequency.

Fig. 7.

Intrinsic modes of single-cell activity of model 1 as a function of NaP andE L. When NaPexceeds 2.2 nS, the model neuron exhibits 3 functional states: silent, oscillatory bursting, and beating. Below 2.2 nS, oscillatory bursting does not occur at any level of E L. Range of NaP used in the figure spans those values used in the simulations presented in this and the companion paper (Butera et al. 1999).

Figure 8 quantifies the effects ofE L on model 2 bursting for a range of values of the subthreshold conductances NaP and KS. For all values of NaP and KS shown, the model possesses regimes of silent, bursting, and beating activity. During bursting activity, the burst period decreases as the neuron is depolarized. However, the parameterization of the dynamics of model 2 byE L differed from that of model 1 in several ways. First, the burst duration increased slightly with depolarization. Second, during tonic firing model 2 could not achieve spike frequencies as low as those obtained in model 1. Third, model 2 supported bursting over a range ofE L (and I app) approximately twice as large as the range of E Lwhere bursting occurred in model 1. These trends persisted as the other parameters of the subthreshold currents (e.g., θk and ςk) were varied. Similar results were obtained when g tonic-e was varied.

Fig. 8.

Mechanisms for frequency control in model 2. Panels show effect of varying E L on burst duration, burst period, burst frequency, and beat frequency of model 2. These curves were calculated for multiple values of NaP (A) and KS (B).I APP below E Lscale shows equivalent variation in I APPwith E L fixed to −65 mV.

The large range of I app in model 2may be attributed to the difference in whole cell conductance (g m) between the two models. Model 2has an additional conductance, KS, which is not present in model 1. This conductance is activated during bursting and beating. The increased g mmeans that the “gain” of the cell is decreased. A larger increment in I app is required to achieve a given increment in depolarization. This is why the frequency-current curves in Fig. 8 are flatter and extend over a largerI app range than those in Fig. 6.

The critical input levels for the transitions between activity modes (i.e., from silence to bursting and bursting to beating), show distinct dependencies on the conductances NaP and KS. In both models (1 and2), the critical value of E L orI app that changes the cell from quiescence to bursting is quite sensitive to NaP. The persistent Na+ current is the primary voltage-dependent current (in model 1, the only one) that begins to activate in the subthreshold regime. At hyperpolarized voltages, these Na+ channels have a low open probability (exponentially small with decreasing V) and the input resistance is high. Thus excitability and the critical level of input needed for bursting is strongly influenced by how many channels there are per unit area, i.e., by NaP. In model 2, KS affects this critical level only modestly, partly because the driving force for KS is much lower than that of NaP and partly because it possesses slow (compared with I NaP) activation kinetics and is relatively inactive at subthreshold voltages. At the bursting-to-beating transition, the cell is depolarized and conductances are well activated. Thus the extra sensitivity to NaP due to high input resistance and low channel activation that we see at the silence to bursting transition is diminished here. The critical levels for E L andI app are almost independent of NaP. The transition is affected by KS because this current controls burst termination (Fig. 5 B3). Although these sensitivities for the silence to bursting transition can be understood qualitatively by considering the subthreshold currents (I sub, see following text) alone, the corresponding sensitivities for the bursting to beating transition are difficult to intuit because the components ofI sub interact with the spike generating currents.

Mechanism for mode transitions and frequency control

In this section, we offer a qualitative biophysical description of how control of the baseline membrane potential, illustrated by varyingE L, alters the activity mode (silence, bursting, and beating as shown in Figs. 4 and 5) of both models. The activity modes exhibited by the model are explained by using steady-state (SS) and quasi-steady-state (QSS) I-V curves of the subthreshold currents (I subI L +I NaP-h for model 1,I subI L+I NaP + I KS formodel 2). This explanation suffices because without the action potential currents I Na andI K, reduced models having only theI sub currents display oscillatory and nonoscillatory states (hyperpolarized silence, oscillatory activity, depolarized silence) that are qualitatively similar to those in the corresponding full models. Figure 9illustrates I-V curves corresponding toI sub in cases shown in Figs. 4 and 5. The SS curves (thick lines) represent I sub versusV when the slow negative feedback process (hinactivation for model 1, k activation formodel 2) is set to steady state [(i.e., h =h (V) in model 1, andk = k (V) inmodel 2]). The QSS curves (thin lines) representI sub versus V when the slow process (h or k) is fixed to a particular value. These curves are generally N-shaped, revealing the regenerative nature of the persistent sodium current. For the remainder of this section, we will only discuss model 1 (Fig. 9 A, 14). An analogous explanation suffices for model 2.

Fig. 9.

Subthreshold currents of model 1(I NaP-h + I L) (left) and model 2(I NaP + I KS +I L) (right) asE L is varied. Values ofE L correspond to values used in Figs. 4 and5 and are −65, −60, −57.5, and −54 mV (A, 1–4) and −65, −59.5, −50, and −40 mV (B, 1–4). Thick lines are steady-state I-V curves [h =h (V) in A, 1–4, k = k (V) inB, 1–4]. Thin I-V plots are quasi-steady-state (QSS) I-V curves, where the slow variable h (A, 1–4) or k(B, 1–4) is fixed to a particular value. InA1 and B1, these values are indicated on the figure. In A, 2 and 3, and B, 2 and 3, these values correspond to the values of h or k at the beginning (a) and end (b) of the active phase of the burst cycle (I andII, top). In A4 and B4,the QSS I-V curve was generated by settingh or k to the mean value during tonic firing (see Figs. 4 B4 and 5B4).

For all values of E L the SS I-Vcurves have a single positive-sloped crossing of the 0 nA axis. This crossing is an equilibrium state that may be stable or unstable. If stable, the equilibrium state corresponds to the model’s resting potential. However, stability is not governed only by the slope’s sign at the zero crossing. If all voltage-dependent conductances vary on time scales faster than the membrane time constant (C M/g m), then a positive slope implies stability. When the time scales are significantly disparate (such as h, which has a time constant of seconds), SS I-V curves may not correctly predict stability, although stability may be inferred through the use of QSS I-V curves. When the QSS I-V curve has a positive crossing of the 0 nA axis, the model is stable at that membrane potential at that particular value of h (i.e., as if h was clamped at the value). Neither the SS I-V nor QSS I-Vcurves provide any information on the dynamics of h itself.

When E L = −65 mV, the model cell is silent (as shown in Fig. 4 A1). The rest potential, at approximately −62 mV where the SS I-V curve ofI sub crosses the 0 nA axis, is stable. The equilibrium point is stable because I NaP-h is relatively deactivated and I NaP-h cannot be adequately recruited to sustain bursting activity. In this case the QSSI-V curves (Fig. 9 A1) show that for all values ofh, even when inactivation is completely removed (h = 1), the hyperpolarized zero crossing is stable (i.e., a positive slope).

For E L equal to −60 mV (Fig. 4 A2) or −57.5 mV (Fig. 4 A3), the model is in the bursting mode and the QSS I-V curves change during the cycle (Gola 1974). An animated QuickTime movie of several burst cycles, illustrating V(t), h(t),I sub(t), and QSSI-V(V, t), is available on the World-Wide Web at Figure9 A, 2 and 3, showsI sub at the maximum (a) and minimum (b) values of h during one complete burst cycle. The burst begins at a value of h (a) whereI sub is inward at all hyperpolarized potentials below the action potential threshold (approximately −45 mV). The cell depolarizes across the action potential threshold and repetitive firing occurs. With each action potential, h inactivates an additional small amount, and this inactivation decreases the inward contribution of I NaP-h toI sub, shifting the QSS I-V curve outward. Firing persists until I sub is net outward at V min, the minimum membrane potential encountered during an action potential. At this point (b), firing stops, and V falls to the hyperpolarized equilibrium point (0 nA crossing with positive slope) of the QSS I-Vcurve (b). As h gradually deinactivates, the QSSI-V curve moves inward, depolarizing the cell as the pseudosteady equilibrium point drifts rightward. Another burst begins when I sub is net inward at all subthreshold potentials (a), and the hyperpolarized equilibrium point of the QSS I-V curve disappears.

As the cell is depolarized further (from A2 toA3, Fig. 4), the values of h encountered during one burst cycle are biased toward smaller values, i.e., less inward current from I NaP-h is necessary to counterbalance the reduced outward current from shiftingE L. In addition, the dynamic range ofh is reduced, from a Δh of ∼0.1 in Fig.4 A2 to <0.02 in Fig. 4 A3. This is the essence of the voltage-dependent frequency control in the model: as the cell is depolarized toward higher values of E L, the cell does not hyperpolarize as far at the end of each burst cycle, and less time is necessary to recover from inactivation ofI NaP-h to begin the following burst. At all values of E L where bursting occurs,I NaP-h is active throughout the burst cycle, and it is the cyclical variation in h that controls the cycle timing.

With still further depolarization such that the equilibrium point of the SS I-V is above the action potential threshold, the model changes to a state of tonic firing (Fig. 4 A4, E L = −54 mV). For this activity mode, h is nearly constant, oscillating at low amplitude with the action potential frequency about a mean value. The QSS I-V curve in Fig.9 A4, generated for this mean value of h (0.315), illustrates that on average I sub is net inward at all subthreshold potentials during tonic spiking.

In summary, the model may be thought of as having two voltage-dependent thresholds—a threshold for burst initiation (approximately −60 mV) and a threshold for action potential initiation (approximately −45 mV). When the equilibrium point of the SS I-V curve is below the burst threshold, the model is quiescent. When the equilibrium point is above the action potential threshold, the model is firing tonically. When the equilibrium point is between these two values, oscillatory bursting is likely. This explanation is only qualitative, however, and a more quantitative explanation may be provided using bifurcation theory (e.g., Rinzel 1985).

Model predictions

In this section, we consider several other functionally important features of the models’ dynamic behaviors that we predict should be observed experimentally.

When action potentials are blocked, each model still exhibits activity modes and transitions that are analogous to those of our full model. Figure 10 demonstrates the models’ dynamic responses when Na = 0. This predicts what we expect in experiments when I Nais blocked (e.g., with QX-314 or TTX–seediscussion), and cells are depolarized progressively. For both models, a subthreshold oscillation remains, with the trends in burst duration, burst period, and interburst depolarization consistent with those of the more complete models (Figs. 4 A, 1–4, and5 A, 1–4). Thus in spite of complex interactions between the spike-generating currents and the subthreshold burst-generating currents, especially at the transition between bursting and beating, these figures demonstrate that the subthreshold currents generally control the burst period and burst duration during repetitive bursting.

Fig. 10.

Oscillatory model behavior in the absence of fast Na+currents. Panels show effects of setting Na to 0 in model 1(A, 1–5) and in model 2 (B, 1–5). Parameter sets for A, 1–4, and B, 1–4, correspond to parameters used in Figs. 4 A, 1–4, and 5 A, 1–4, respectively. Parameters for A5 and B5 areE L = −53 mV andE L = −32 mV, respectively. - - -, −60 and −30 mV references.

Cells in hyperpolarized silent mode may retain burst excitability if NaP is not too small, showing single burst responses to transient inputs. A brief depolarizing input (50 ms) of sufficient magnitude is capable of triggering a single sustained burst lasting several hundred milliseconds (Fig.11 A). Similar responses occur across a wide range of values of E L and NaP where the model is initially silent, providing that NaP is not too small. At rest, I NaP-h is relatively deactivated (m is low) and deinactivated (h is high). A brief depolarization activates I NaP-h. BecauseI NaP-h inactivates slowly,I NaP-h remains relatively deinactivated and a larger conductance for I NaP-h is elicited, generating sufficient inward current to trigger a single burst. As analyzed in the companion paper (Butera et al. 1999), such burst triggering is important for recruiting inactive cells and achieving synchronous bursting in a synaptically coupled population of these pacemaker cells with heterogeneous properties where a substantial fraction of the population is intrinsically in a silent mode.

Fig. 11.

Transient burst responses. A: brief but sufficiently large transient depolarizing input can elicit a single sustained burst. Model was at rest (E L = −65 mV), and depolarizing current pulses (50 ms, 10 and 15 pA) were applied.B: transient hyperpolarizing input of sufficient magnitude and duration may elicit a single posthyperpolarization rebound burst. Model was at rest (E L = −62 mV), and a hyperpolarizing current pulse (500 ms, 60 pA) was applied. In both cases, the parameters E L and NaP correspond to points in Fig. 7 to the left of the bursting regime.

A single burst also may be elicited as a form of posthyperpolarization rebound (Fig. 11 B). A sustained hyperpolarizing input removes some inactivation from I NaP-h (i.e.,h increases). On abrupt release of the hyperpolarizing input, the model depolarizes. Because inactivation ofI NaP-h develops slowly,I NaP-h has an increased conductance asV nears its rest value, and this increased conductance may trigger a single posthyperpolarization burst. This response may occur across a wide range of parameters provided that the stimulus is of sufficient duration (the inactivation ofI NaP-h is a slow process), the hyperpolarizing stimulus is of sufficient magnitude (so thath increases sufficiently), and the resting neuron is not too hyperpolarized (so that sufficient h remains to be deinactivated). For example, the response of Fig. 11 B was generated for E L = −62 mV. WhenE L = −65 mV, the resting value of his sufficiently large (0.92) that even a large hyperpolarization (to bring h nearly to 1) does not produce an increase in the conductance of I NaP-h strong enough to elicit a posthyperpolarization burst, unless an additional depolarizing input is coincident with the release from hyperpolarization. Thus we predict that this response can be elicited by hyperpolarization release only for a window of membrane potentials that is not too hyperpolarized.

During repetitive bursting, transient inputs can reset the rhythm. A brief (50 ms) hyperpolarizing input (10 pA) is capable of resetting the active phase of a burst (Fig. 12). This resetting reduces the duration of the current burst as well as the following interburst interval. The duration of the following interburst interval varies with the duration of the preceding burst. The duration of the subsequent burst is unaffected. For example, a transient current pulse applied early in the burst (b 1) terminates the burst, and the following interburst interval is significantly shortened, with the next burst firing at b 2. A transient current pulse applied late in the burst (c 1) also terminates the burst, with the next burst firing at c 2. Although the poststimulus interburst interval is shortened in both B and C,the poststimulus interburst interval in B is shorter. This behavior would be important for controlling burst cycle timing under conditions where there is dynamic resetting by inhibitory synaptic mechanisms (see discussions in Smith 1997; Smith et al. 1995).

Fig. 12.

Resetting response. Brief but sufficiently strong transient hyperpolarizing input can reset a burst. Simulations performed withE L = −59 mV, which corresponds to a burst period of ∼4 s. Control period shown in A, with bursts occurring at times α1 and α2. Simulations in B and C were started with identical initial conditions, and a 50-ms, 10-pA hyperpolarizing input was applied early (b1) or late (c1) in the burst cycle. In both cases, the brief input reset the burst with the following burst occurring earlier. Duration until the subsequent burst was proportional to how early or late in the burst the transient input was applied.


Justification of models

We have developed minimal models for oscillatory bursting neurons found in the pre-Bötzinger complex in vitro. Although there is strong evidence that neurons with bursting pacemaker-like properties are generating the respiratory rhythm in vitro, it has not been proven that the neurons on which our models are based, with the behavior depicted in Fig. 1, actually generate the rhythm. Such causality remains difficult to establish definitively experimentally. However, these cells are the only inspiratory neurons with intrinsic pacemaker-like properties that have been identified in the pre-Bötzinger complex from electrophysiological recording and mapping of neuron activity in slice preparations (Johnson et al. 1994; Smith et al. 1991; Koshiya and Smith 1998, 1999). They are currently candidates for the rhythm-generating cells in vitro (see Smith et al. 1995).

Intrinsic conductances responsible for the oscillatory bursting behavior of these cells have not yet been identified experimentally. Furthermore there is likely to be considerable heterogeneity in cellular parameters for different neurons in the rhythm-generating cell population, e.g., maximal conductances of burst-generating currents ( NaP) and baseline membrane potentials (E L); this heterogeneity is treated in the companion paper (Butera et al. 1999). Figure 1 only provides one example of voltage-dependent behavior in a bursting inspiratory neuron. Variations in properties such as the voltage range for oscillatory bursting, oscillatory frequency range, and the extent to which a given neuron in the population will exhibit oscillatory bursting (see Functional states of pacemaker neurons), remain to be quantified experimentally. Given the limited experimental data, we therefore have postulated a minimum set of conductance mechanisms and explored a range of parameter values to determine if they can account for the observed behaviors as well as to investigate the range of behaviors likely to be exhibited by cells with these types of conductances. Besides fast Na+ and K+currents responsible for the generation of action potentials, ionic mechanisms must be postulated for both the initiation and termination of bursting. Although we have assigned specific types of currents to these functions, as justified in the following text, the models are general and are designed to provide insights into plausible mechanisms.


In numerous invertebrate and mammalian bursting neurons, the onset of bursting is caused by an inward cationic current that activates at subthreshold potentials. This current is responsible for maintaining the negative-slope region of the I-V curve. T-type and persistent Ca2+ currents have been implicated in burst generation in mammalian neurons (e.g., Llinás 1988). Several types of Ca2+ currents exist in respiratory neurons (Onimaru et al. 1996); however, inspiratory neurons have been identified which continue to undergo rhythmic bursting under low Ca2+ conditions (Johnson et al. 1994) (see also Fig. 1). This evidence suggests that Ca2+ currents are not essential for the generation of bursting behavior. Therefore although it is known that the pre-BötC oscillatory bursting cells have Ca2+currents (Koshiya and Smith 1998) and it is likely that Ca2+ currents play a role in modifying a neuron’s firing properties, they are not considered in our formulation of a minimal model for bursting activity.

Calcium-activated nonspecific cation (CAN) currents have been implicated in playing a role in oscillatory bursting neurons inHelix (Swandulla and Lux 1985) andAplysia (Kramer and Zucker 1985) and are ubiquitous in a variety of neuronal preparations (see Partridge and Swandulla 1988 for review). Their role has been postulated to maintain the active phase of a burst. Because oscillatory pre-BötC neurons continue to burst in a low Ca2+medium, we do not consider CAN currents necessary to maintain oscillatory bursting. In addition, because CAN currents have conductances that are typically voltage independent, it is unlikely that they could be responsible for generating the rapid depolarization associated with the initiation of a burst.

The mixed cationic current I H also has been implicated in rhythmic bursting activity in a variety of neuronal preparations (e.g., Calabrese and DeSchutter 1992;Lüthi and McCormick 1998; McCormick and Huguenard 1992). However, an I H current is unlikely to be necessary because rhythmic bursting may occur even when there is very little hyperpolarization during the interburst interval (Fig. 1). Although I H may not be essential for bursting, this current could be involved in regulating the time course of membrane depolarization during the interburst interval and therefore would be important for regulating burst frequency. A relatively flat membrane potential trajectory could occur in the presence of I H if this current was balanced by an outward K+-dominated leakage current, as incorporated in our model, or another K+ leak current such as an inward rectifier (I R).I R has been shown in other preparations to be a modulatory target that regulates the frequency of bursting (Benson and Levitan 1983; Butera et al. 1995; Drummond et al. 1980) and has been proposed to be involved in the frequency control of pre-BötC pacemaker cells (Smith et al. 1995). More experimental information is required on the mix of subthreshold currents controlling the time course of the interburst interval. Our model was designed to investigate a minimal set of essential burst-generating conductance mechanisms.

In light of these observations, a compelling alternative remains for burst initiation: inward cationic currents that activate at subthreshold membrane potentials. For example, a slow-inward Ca2+ current has been reported in Aplysia neuron R15 that has voltage-dependent activation and a slower Ca2+-dependent inactivation (Adams and Levitan 1985). Because bursting persists under low-Ca2+conditions, another plausible candidate is the persistent Na+ current (I NaP).I NaP has been implicated in bursting behavior in rat hippocampal, hypothalamic, and cortical neurons (Franceschetti et al. 1995; Li and Hatton 1996; Llinás 1988) and identified in a variety of other mammalian neurons. The voltage dependence ofI NaP has been quantitatively characterized in mouse and rat neocortical neurons (Brown et al. 1994;Fleidervish et al. 1996; French et al. 1990). We adopted I NaP as the current responsible for burst initiation in our model. Ongoing experiments in our laboratory seek to identify the existence ofI NaP or a similar subthreshold cationic current in pre-BötC inspiratory bursting neurons.

French et al. (1990) and Brown et al. (1994) report values for the parameters of the voltage-dependent activation of I NaP: θm = −50 and −48 mV, respectively, and absolute values of ςm in the range of 4–9 mV. Our choice of θm (−40 mV) is depolarized from these values. However, similar voltage-dependent frequency control (at more hyperpolarized potentials) is obtained from the model if θm is reduced to −48 mV with a similar shift in θh . The formulation of m (V) forI NaP in our model is consistent with observations from cortical neurons where the threshold forI NaP activation is ≥10 mV hyperpolarized to the threshold for I Na activation (Brown et al. 1994; Stafstrom et al. 1982) and is present in a voltage range where most other voltage-gated currents are inactive. Thus I NaP dominates the excitability of the membrane at subthreshold potentials (Stafstrom et al. 1982). This is evident when comparingm(V) for I Na withm (V) forI NaP in Fig. 2 (A and B1).


Spike-frequency histograms of pre-BötC bursting neurons in rhythmically active slices (Johnson et al. 1994) show that the firing frequency decays throughout the burst phase (Fig. 1). Also the initiation and termination of the burst are accompanied by a rapid transition between the silent phase and the firing of action potentials (Fig. 1) and vice versa. In light of theoretical studies (Bertram et al. 1995a; Rinzel 1987) of the mechanisms by which individual neurons may exhibit bursting behavior, both of these lines of evidence suggest that a minimal mechanism for bursting requires just one slow recovery process, such as intracellular Ca2+-accumulation or a slow voltage-dependent gating mechanism. Because rhythmic bursting pre-BötC neurons have been identified in low-Ca2+ preparations (Johnson et al. 1994), we postulate that the recovery process is voltage dependent. This process could be either a slowly inactivating inward current or a slowly activating outward current (presumably K+).

Few inward cationic currents have been identified that have voltage-dependent kinetics on the order of seconds that are required to produce the bursting dynamics of the pre-BötC neurons. Slow voltage-dependent inactivation of a Ca2+-current has been identified in pancreatic β-cells with a time constant of 2–6 s (Hopkins et al. 1991), but a similar current has not been identified in neuronal preparations. A persistent Na+current is the main candidate, but the slow inactivation parameters ofI NaP-h have not been as well characterized; the parameters used in our models are consistent with available measurements. Our h (V) is similar in shape to the voltage dependency of the slow inactivation ofI NaP-h reported by Fleidervish et al. (1996). The steady-state I-V plots ofI NaP-h in our model (Fig. 2 B2) are consistent in shape with the I-V characteristics estimated by Fleidervish and Gutnick (1996) during slow (2.33 mV/s) and fast (70 mV/s) voltage ramps. When using ramps of ≥35 mV/s, they reported that I NaP begins to activate around −60 mV and reaches a peak by −25 mV. Fleidervish and Gutnick (1996) reported the time constant for the onset of slow inactivation of I NaP-h as 2.06 s at +20 mV and the time constant for recovery from slow inactivation ofI NaP-h as 2.31 s at −70 mV and 1.10 s at −90 mV. These points are within the range of our model of τh(V) (Fig. 2 B1), although they lie at extreme ends of the bell-shaped curve τh(V). We used a value of 10 s forτ̄h to produce burst periods on the time scales that are observed experimentally.

The hypothetical burst termination mechanism in model 2requires a slowly activating voltage-dependent K+ current. Slow voltage-gated Phe-Met-Arg-Phe-NH2 (FMRF-amide) sensitive K+-currents (I KF) recently have been identified in heart interneurons in the leech (Nadim and Calabrese 1997). This current has a time constant of several seconds. However, K+ currents with similar voltage-dependent kinetics have not yet been identified in mammalian neurons.

Although we have postulated voltage-dependent slow processes for burst termination, it is possible that these processes may be Ca2+ dependent as well because a Ca2+ gradient still may exist under the low-Ca2+ conditions where bursting behavior has been observed (Johnson et al. 1994). However, our experiences with the dynamic mechanisms of these classes of models (see Comparison with other bursting cells) is that dynamically similar mechanisms for bursting may be obtained whether the slow process is voltage or Ca2+dependent.

Comparison of models

Models 1 and 2 both demonstrate several qualitative features of the oscillatory bursting behavior of pre-BötC inspiratory neurons. These features include voltage-dependent control of the model activity: as the neuron is depolarized via I app org tonic-e, or varying E L(e.g., by controlling [K+]o), the neuron can be biased from silence to bursting to beating activity. During bursting, the spike-frequency decays during the burst and the bursting frequency is strongly dependent on the level of depolarization. However, model 1 is more consistent with our experimental data in two respects: the membrane potential trajectory during the interburst interval is much flatter (in membrane potential) and the burst duration decreases as the baseline membrane potential is biased toward more depolarized levels. It is possible that a more complex repertoire of ionic currents may allow model 2 to possess a flatter interburst interval (e.g., see Smith et al. 1995). Also, it is possible that we have not identified a possible parameter set for model 2 where the burst duration decreases, rather than increases, with depolarization. Nevertheless, we have searched parameter space for all of the parameters of the subthreshold currents (specifically, maximal conductances and parameters responsible for the dynamics of the gating variables: KS, Nap, θm, ςm, θk, ςk) in model 2 and have not identified parameter regimes wheremodel 2 exhibits such behavior.

We found that by varying additional parameters ofI KS we can mathematically transform the ionic current equations of model 1 into a “nearly dynamically equivalent” model that possess the ionic current equations ofmodel 2. Although the nature of this transformation is beyond the scope of this paper, the nearly equivalent model 2 possesses a relatively flat interburst interval, and a burst duration that decreases with depolarization until near the threshold for beating. However, to achieve these dynamics using the ionic currents of model 2, we found it necessary to either make the reversal potential of I KS −55 mV or incorporate an additional fast-activating gating variable toI KS such that I KS is only active at membrane potentials above −55 mV. With either of these manipulations to model 2, g mdecreases throughout almost the entire silent phase.

The relative flatness of model 1’s interburst voltage time course can be understood as follows. During the interburst phase, the subthreshold currents are essentially balanced and they sum to zero:I sub ≈ 0. Thus V is essentially at a hyperpolarized pseudo-steady-state and drifts upwards as the gating variable [h(t) or k(t)] slowly evolves. By differentiating this pseudo-steady-state relation, one can estimate the expected variations in V as the slow gating variable evolves: ΔVh or ΔVk. For model 1, we find that ΔVh is proportional toI NaP (i.e., I NaP-h withh = 1), while for model 2ΔVk is proportional to KS(VE K) (i.e., I KS withk = 1). Because I NaP-h is almost totally deactivated (m ≈ 0) in the voltage range of the interburst phase we see that ΔVh is much less than ΔVk. This qualitative prediction of model 1’s flat interburst voltage trajectory was partial motivation for us to develop this mechanism early in our study. A similar argument shows that the interburst voltage is more sensitive to I app in model 1 than inmodel 2 by a factor of ≥2.

Models with a voltage- and time-dependent burst inactivation process will exhibit some depolarizing trajectory of the interburst membrane potential, as shown by model 1 and to a greater extent bymodel 2. It has been speculated that postburst hyperpolarization and depolarizing drift of interburst membrane potential may be a signature of pre-BötC pacemaker cells (seeRekling and Feldman 1998). Recordings of the candidate pre-BötC pacemaker neurons, however, can exhibit nearly flat interburst interval trajectories (Fig. 1) (Smith et al. 1991). This may indicate the presence of a more complex mix of subthreshold currents (e.g., Smith et al. 1995) and/or a nonuniform spatial distribution of ion channels between the soma and proximal dendrites (e.g., Li et al. 1996).

Experiments to test models

In light of the differences discussed above and the robustness of the dynamics of model 1, we are inclined to presentmodel 1 as a more feasible mechanism for oscillatory bursting in pre-BötC neurons. However, due to the limitations previously discussed, we must test for both mechanisms in vitro. One test that could support or refute the ionic currents of model 1 would be to measure the membrane conductance at various points throughout the silent phase of the cycle (Atwater and Rinzel 1986). In model 1, the membrane conductance increases monotonically as I NaP-h recovers from inactivation (Fig. 4 B, 1–4). The recovery of hand the gradual depolarization of the membrane potential (via the fast activation of I NaP-h) both contribute to an increase in Nap-h, and thusg m . The presence of additional subthreshold currents with voltage-dependent conductances (e.g.,I R or I H), however, could alter this trend.

In model 2 during the silent phase of the burst cycle, there are two competing effects with respect tog m . The deactivation ofI KS contributes toward a reduction ofg m , while the membrane’s gradual depolarization activates I NaP, contributing toward an increase in g m . Thus the overall change in g m in model 2depends on the relative strength of I NaP andI KS. In Fig. 5 it is shown thatg m decreases then increases during the silent phase of the burst cycle.

Another test to differentiate between the two mechanisms would be to apply a voltage-clamp prepulse to 0 mV for 10 s and then hyperpolarize to −75 mV. In model 2, the prepulse would activate I KS, and on hyperpolarization (as long as we hyperpolarize above E K) we would expect to see a slowly decaying outward current. In model 1 we would not expect to see the slow decay as long as the clamp potential was below the activation threshold for I NaP-h. Therefore the existence of a slow current following a long depolarized prepulse would demonstrate sufficiency for model 2 and would be inconclusive for model 1.

Figure 6 B quantifies the effect of a mean level of tonic synaptic drive on the behavior of model 1. This cannot be tested directly in vitro, since the act of decoupling excitatory synaptic connections between bursting neurons (e.g., via a low Ca2+ bathing medium) also blocks tonic synaptic input to these cells as well. However, this mechanism could be demonstrated through the use of a dynamic clamp (Sharp et al. 1993) after synaptic connections have been blocked. Nevertheless these results illustrate that the depolarizing inputs of extrinsic tonic drive have similar effects on the activity of model 1 as varying E L or I app. Similar results to Fig. 6 are obtained as E K is varied (with I L separated into Na+and K+ components). We have found that [K+]o regulates the oscillation frequency of individual pacemaker-like cells under low-Ca2+ or 6-cyano-7-nitroquinoxaline-2,3-dione conditions (R. Butera, C. Del Negro, and J. Smith, unpublished data). Similar results to those shown in Fig. 8 for model 2 are obtained asE K or g tonic-e is varied.

The subthreshold mechanisms responsible for oscillatory bursting in pre-BötC neurons may be elucidated more easily experimentally if the fast Na+ current is blocked, as in Fig. 10. Our models predict that a subthreshold oscillation should remain under these conditions, allowing differences in the trajectories and voltage dependence of potential to be distinguished for the two models.I NaP, in some mammalian neurons, is TTX-sensitive (Fleidervish and Gutnick 1996), necessitating block of fast Na+ by intracellular QX-314, although TTX-insensitive forms of I NaP also have been described (Hoehn et al. 1993; Oka 1996).

It also will be necessary for recordings from oscillatory bursting pre-BötC neurons to possess transient responses similar to those shown in Figs. 11 and 12. These tests will not provide evidence in favor of model 1 or model 2 but will simply demonstrate that either model is consistent with experimental data. The responses are not unique, however, and a more complex bursting neuron with additional slow variables and ionic currents may possess similar responses (e.g., Butera et al. 1995, 1997a; Demir et al. 1997). Rebound bursting has been proposed as an important mechanism for producing transitions from inactive to active phases in neural oscillators, particularly when there is phasic inhibition hyperpolarizing burst-generating cells. I H andI T have been identified as ionic mechanisms that contribute to rebound bursting, and these currents may exist in respiratory neurons (Ramirez and Richter 1996;Richter 1996). Our models demonstrate thatI NaP-h is another current mechanism promoting posthyperpolarization bursting (Fig. 11 B). This may be functionally important in the respiratory oscillator under conditions where the pacemaker cells are embedded in the respiratory pattern generation network and receive phasic inhibitory inputs (see discussions in Smith 1997).

Comparison with other bursting cells

Voltage-dependent frequency control of oscillatory bursting has been described in a variety of neuronal preparations, including neuron R15 in Aplysia (Mathieu and Roberge 1971;Wilson 1982), neuron AB in stomatogastric ganglion (Abbott et al. 1991), medial mammalian body (MMB) neurons in the guinea pig hypothalamus (Alonso and Llinás 1992), and magnocellular neurons in the rat hypothalamus (Li and Hatton 1996). In all of these cases, the frequency of bursting increases, and the baseline membrane potential is depolarized as I app is increased. Similar effects of I app on burst frequency are exhibited by models of R15 (Canavier et al. 1991) and AB (Abbott et al. 1991; Epstein and Marder 1990).

The models presented in this paper often are described as square-wave, or type I, bursters (Bertram et al. 1995a; Rinzel 1987). Various models of this type have been used to describe the electrical bursting activity of pancreatic β-cells (Bertram et al. 1995b; Chay and Keizer 1985; Keizer and Smolen 1991; Sherman and Rinzel 1992). In all of these models, burst termination occurs by either a slowly activating K+ current or a slowly inactivating Ca2+ current. Models of electrical bursting in pancreatic-β cells that rely on either of these slow processes have used time constants with maximal values of 30–50 s (Bertram et al. 1995a; Keizer and Smolen 1991;Sherman and Rinzel 1992). The effect of membrane depolarization in these secretory cells has been studied by modifying a glucose-dependent conductance and not with I app. Thus we reexamined a few of these β-cell models (Sherman and Rinzel 1992; Sherman et al. 1990; Smolen and Keizer 1992) to study the effects ofI app and found that all three demonstrate voltage-dependent frequency control as described in the preceding text. To our knowledge, the models that we developed here are among the first examples of type I bursting (Bertram et al. 1995a;Rinzel and Lee 1987) in a neuronal preparation. The appearance of type-I like bursting behavior has also been reported in trigeminal motoneurons (Del Negro et al. 1998).

It is difficult to generalize any principles regarding voltage-dependent control of burst duration. In some cases, the burst duration increases with I app, such as R15 (Canavier et al. 1991; Mathieu and Roberge 1971; Wilson 1982) and MMB neurons (Alonso and Llinás 1992), whereas in other cases burst duration decreases (Sherman and Rinzel 1992;Sherman et al. 1990) or has very little change (Abbott et al. 1991; Epstein and Marder 1990). Models with similar general mechanisms for burst initiation and termination (e.g., model 2 of this paper) (Sherman and Rinzel 1992) yield opposite results with respect to the voltage-dependent control of burst duration. These results may depend on the interaction of action potential currents with subthreshold processes, a dynamically complex process (de Vries and Miura 1998; Pernarowski et al. 1992;Terman 1991) that warrants further theoretical investigation.

Intrinsic control mechanisms

Neuromodulatory afferent inputs to the pre-BötC pacemaker cells that regulate maximum conductances or voltage-dependent parameters of the intrinsic subthreshold currents would play a major role in burst frequency/duration control (see discussions inSmith 1997; Smith et al. 1995).Models 1 and 2 offer uniquely different possibilities with respect to the possible role of subthreshold currents as neuromodulatory targets. For model 1ḡ NaP may be viewed as a tunable parameter that scales the dynamic range of bursting (Figs. 6 and 7). As NaP is increased, two effects are observed: the range of values of E L where bursting is supported is increased and the range of burst periods possible as E L is varied is increased. The burst duration is primarily a function of EL and relatively independent of values of NaP.

Model 2 offers two distinct subthreshold conductances for controlling the model’s burst dynamics, NaP and KS(Fig. 8). Decreasing KS and/or increasing NaP increases burst duration. KS and NaPhave complementary roles with respect to controlling the frequency of bursting: NaP controls the minimal value of E L where bursting is supported and the maximum burst period, whereas KS controls the maximal value of E L (org tonic-e or I app) where bursting is supported and the minimal burst period. Although KS appears to set an upper bound on the range of E L values that support bursting, Figure8 B2 illustrates that KS has a minimal effect on the overall burst period, which may be considered to be primarily a function of E L. Thus KS may serve as a mechanism for varying duty cycle with minimal effects on the overall burst frequency.

In both models, g L provides a mechanism for controlling baseline membrane potential and burst frequency. K+-dominated conductances that provide a leakage or background current in the subthreshold voltage range, including inward rectifying K+ conductances, have been proposed to play an important role in control of pre-BötC pacemaker cell activity (Smith et al. 1995) as a target for neuromodulation. Our models show that high g L results in quiescence. As g L is reduced, bursting occurs at a frequency dependent on the value of g L. Sufficiently low values of g L cause beating. Thus the K+-dominated leak conductance can control the activity mode (quiescence, bursting, beating) and provide frequency control in the voltage range for bursting.

Functional states of pacemaker neurons

The model pacemaker neurons have several functional states as a consequence of the voltage-dependent properties ofI NaP: quiescence, oscillatory bursting, and beating. Previously we have called these pacemaker cells conditional pacemakers to signify that conditions (appropriate ranges for level of depolarization or magnitude of burst-generating conductances) must exist for oscillatory bursting to occur (Smith et al. 1991,1995). As indicated in Figs. 6 and 7 for model 1, the magnitude of NaP determines whether the neuron will exhibit oscillatory bursting: at low values of NaP the cells only exhibit quiescence or beating (Fig. 7), but they still (but not necessarily) may be burst capable so that a self-terminating burst can be triggered by a brief transient depolarization. As NaP is increased, the cells exhibit a voltage regime where intrinsic oscillatory bursting occurs. In our pacemaker cell population models, as analyzed in detail in the companion paper (Butera et al. 1999), we presume that there is ordinarily heterogeneity in values of NaP andE L within the cell population such that there is a distribution of cells in the different states; only a fraction of the cells (when uncoupled) are in the oscillatory bursting mode. With excitatory synaptic coupling, synchronous oscillatory bursting occurs and this heterogeneity results in a functionally more robust oscillator than if the cells were homogeneous (e.g., all in the oscillatory bursting mode). Because of the inactivation properties ofI NaP and synaptic interactions, synchronous oscillations can occur in the coupled population even if none (lowI NaP) or few of the cells exhibit oscillatory bursting. As discussed earlier, NaP andg L are parameters potentially tunable by neuromodulation. We therefore view the range of behaviors exhibited by the pacemaker cells as a functional continuum, and the fraction of cells in the population that are in any state could be tunable by neuromodulatory afferent inputs.


The minimal models we have developed provide plausible mechanisms for generating the multistate, voltage-dependent behavior observed for the candidate pre-BötC pacemaker neurons. We conclude that a burst-generating mechanism dominated by the activation-inactivation properties of a single cationic conductance, postulated to be a persistent Na+ conductance, can account for the voltage-dependent bursting behavior. The dynamics of the burst cycle are controlled principally by the kinetics of inactivation and recovery from inactivation of this conductance. Although the actual burst-generating current in the pre-BötC pacemaker neurons remains to be identified, our model indicates the essential features that are required to produce the experimentally observed bursting behavior.


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

This work was supported by the intramural research programs of the National Institutes of Health.

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


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

  • The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.


View Abstract