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

Robert J. Butera, Jr.,1,2 John Rinzel,1,2,3 and Jeffrey C. Smith1

 1Cellular and Systems Neurobiology Section, Laboratory of Neural Control, National Institute of Neurological Disorders and Stroke, National Institutes of Health;  2Mathematical Research Branch, National Institute of Diabetes and Digestive and Kidney Diseases, National Institutes of Health, Bethesda, Maryland;  3Center for Neural Science and Courant Institute of Mathematical Sciences, New York University, New York City, New York 10013


    ABSTRACT
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

Butera, Robert J., Jr., John Rinzel, and Jeffrey C. Smith. Models of Respiratory Rhythm Generation in the Pre-Bötzinger Complex. I. Bursting Pacemaker Neurons. J. Neurophysiol. 82: 382-397, 1999. 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 INaP-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.


    INTRODUCTION
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

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


    METHODS
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

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 http://netlib.cs.utk.edu/ode/cvode.tar.Z. 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 ftp://ftp.math.pitt.edu/pub/bardware.

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 also Koshiya 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. 1B), 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. 1B) (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.



View larger version (30K):
[in this window]
[in a new window]
 
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 in A 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 (INaP). However, burst termination in the two models occurs by contrasting mechanisms. In model 1, burst termination occurs by the slow-inactivation of INaP-h, a persistent Na+-current with slow inactivation. In model 2, a slowly activating K+-current (IKS) 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.

FORMULATION OF MODEL 1. 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:
<IT>C</IT> <FR><NU><IT>d</IT><IT>V</IT></NU><DE><IT>d</IT><IT>t</IT></DE></FR><IT>=</IT>−<IT>I</IT><SUB><IT>NaP</IT></SUB><IT>−</IT><IT>I</IT><SUB><IT>Na</IT></SUB><IT>−</IT><IT>I</IT><SUB><IT>K</IT></SUB><IT>−</IT><IT>I</IT><SUB><IT>L</IT></SUB><IT>−</IT><IT>I</IT><SUB><IT>tonic-e</IT></SUB><IT>+</IT><IT>I</IT><SUB><IT>app</IT></SUB> (1)
where 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 to
<FR><NU>d<IT>x</IT></NU><DE><IT>d</IT><IT>t</IT></DE></FR><IT>=</IT><FR><NU><IT>x</IT><SUB><IT>∞</IT></SUB>(<IT>V</IT>)<IT>−</IT><IT>x</IT></NU><DE><IT>&tgr;</IT><SUB>x</SUB>(<IT>V</IT>)</DE></FR> (2)

<IT>x</IT><SUB><IT>∞</IT></SUB>(<IT>V</IT>)<IT>=</IT>{<IT>1+exp</IT>[(<IT>V</IT><IT>−&thgr;</IT><SUB>x</SUB>)<IT>/&sfgr;</IT><SUB>x</SUB>]}<SUP><IT>−1</IT></SUP> (3)

&tgr;<SUB>x</SUB><IT>=<A><AC>&tgr;</AC><AC>&cjs1171;</AC></A></IT><SUB>x</SUB><IT>/cosh </IT>[(<IT>V</IT><IT>−&thgr;</IT><SUB>x</SUB>)<IT>/</IT>(<IT>2&sfgr;</IT><SUB>x</SUB>)] (4)
where xinfinity (V) is the steady-state voltage-dependent (in)activation function of x and tau x(V) is the voltage-dependent time constant. xinfinity (V) is a sigmoid with a half-(in)activation at V = theta x and a slope that is proportional to 1/sigma x(V) is a bell-shaped curve that has a maximal value <A><AC>&tgr;</AC><AC>&cjs1171;</AC></A>x at V = theta x and a half-width determined by sigma x. Thus each gating variable is described by only three parameters. The preceding formulation also may be expressed in the HH formalism dx/dt = alpha x(V)(1 - x- beta x(V)x, where alpha x(V) = [1/(2<A><AC>&tgr;</AC><AC>&cjs1171;</AC></A>x)]exp[-(V - theta x)/2sigma x)] and beta x(V) = [1/(2<A><AC>&tgr;</AC><AC>&cjs1171;</AC></A>x)]exp[(V - theta x)/(2sigma x)]. This is the simplest possible formulation for the voltage dependency of alpha (V) and beta (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 (INa) and a delayed-rectifier K+ current (IK). The equations for these currents are inspired by the HH-formulation, but the gating variables satisfy a reduced voltage-dependent description
<IT>I</IT><SUB><IT>Na</IT></SUB><IT>=</IT><IT><A><AC>g</AC><AC>&cjs1171;</AC></A></IT><SUB><IT>Na</IT></SUB><IT>m</IT><SUP><IT>3</IT></SUP><SUB><IT>∞</IT></SUB>(<IT>V</IT>)(<IT>1−</IT><IT>n</IT>)(<IT>V</IT><IT>−</IT><IT>E</IT><SUB><IT>Na</IT></SUB>) (5)

<IT>I</IT><SUB><IT>K</IT></SUB><IT>=</IT><IT><A><AC>g</AC><AC>&cjs1171;</AC></A></IT><SUB><IT>K</IT></SUB><IT>n</IT><SUP><IT>4</IT></SUP>(<IT>V</IT><IT>−</IT><IT>E</IT><SUB><IT>K</IT></SUB>) (6)
where the parameters of INa are <A><AC>g</AC><AC>&cjs1171;</AC></A>Na = 28 nS, ENa = 50 mV, theta m = -34 mV, and sigma m = -5 mV and the parameters of IK are <A><AC>g</AC><AC>&cjs1171;</AC></A>K = 11.2 nS, EK = -85 mV, theta n = -29 mV, sigma n = -4 mV, <A><AC>&tgr;</AC><AC>&cjs1171;</AC></A>n = 10 ms. The formulation for INa 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 as n and is approximated by h = 1 - n (Krinskii and Kokoz 1973; Rinzel 1985). The voltage-dependent functions for minfinity 3(V) and hinfinity (V) for INa and ninfinity 4(V) and tau n(V) for IK are illustrated in Fig. 2A. Quantitative measurements of the gating characteristics of INa and IK 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 for INaP-h in model 1 is identical to that for INaP in model 2 except that INaP-h also has the inactivation term h. The following discussion regarding the activation parameters of INaP-h applies to INaP (in model 2) as well. INaP-h is of the form
<IT>I</IT><SUB><IT>NaP-h</IT></SUB><IT>=</IT><IT><A><AC>g</AC><AC>&cjs1171;</AC></A></IT><SUB><IT>NaP</IT></SUB><IT>m</IT><SUB><IT>∞</IT></SUB>(<IT>V</IT>)<IT>h</IT>(<IT>V</IT><IT>−</IT><IT>E</IT><SUB><IT>Na</IT></SUB>) (7)
where <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP = 2.8 nS, theta m = -40 mV, sigma m = -6 mV, theta h = -48 mV, sigma h = 6 mV, <A><AC>&tgr;</AC><AC>&cjs1171;</AC></A>h =10,000 ms, and ENa = 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) of INaP-h is distinct from that of INaP, where similar variable names were used. Throughout the rest of this text, any reference to h will always refer to the inactivation of INaP-h and not INa.



View larger version (23K):
[in this window]
[in a new window]
 
Fig. 2. Gating and I-V characteristics of components of models 1 and 2. A: spike generating kinetics. minfinity 3(V) and hinfinity (V) of INa and ninfinity (V) and tau n(V) of IK: note that h = 1 - n. B1: gating characteristics of INaP: minfinity (V), hinfinity (V), and tau h(V) (bold). Left: y axis scale for steady-state gating functions. Right: y axis scale for tau h(V). B2: I-V plots of INaP for h = hinfinity (V) and h = 1. First case results in a small window current at subthreshold potentials. Second case corresponds to INaP-h with complete removal of inactivation. C1: gating characteristics of IKS: kinfinity (V) and tau k(V) (bold). Left: y axis scale for activation function; right: y axis scale for tau k(V). C2: I-V plots of INaP + IKS for k = kinfinity (V) and k = 0. First case results in a small current at subthreshold potentials. Second case corresponds to INaP with complete removal of the opposing IKS.

Figure 2B1 illustrates the functions minfinity (V), hinfinity (V), and tau h(V) for INaP-h. Figure 2B2 illustrates INaP-h when h hinfinity (V) and h = 1. This figure predicts the INaP-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 IL, defined as IL = <A><AC>g</AC><AC>&cjs1171;</AC></A>L(V - EL). This current is K+ dominated. The nominal parameters for IL are <A><AC>g</AC><AC>&cjs1171;</AC></A>L = 2.8 nS and EL=-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 conductance gtonic-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 Itonic-e gtonic-e(V - Esyn-e), where Esyn-e is 0 mV, the reversal potential of non-NMDA glutamatergic synaptic currents. Throughout this paper Itonic-e is considered to be inactive (gtonic-e = 0 nS) unless otherwise indicated. Itonic-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 Iapp is included for completeness and is set to 0 pA unless otherwise noted.

FORMULATION OF MODEL 2. 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, IKS. Thus the current balance equation is
<IT>C</IT> <FR><NU><IT>d</IT><IT>V</IT></NU><DE><IT>d</IT><IT>t</IT></DE></FR><IT>=</IT>−<IT>I</IT><SUB><IT>NaP</IT></SUB><IT>−</IT><IT>I</IT><SUB><IT>KS</IT></SUB><IT>−</IT><IT>I</IT><SUB><IT>Na</IT></SUB><IT>−</IT><IT>I</IT><SUB><IT>K</IT></SUB><IT>−</IT><IT>I</IT><SUB><IT>L</IT></SUB><IT>−</IT><IT>I</IT><SUB><IT>tonic-e</IT></SUB><IT>+</IT><IT>I</IT><SUB><IT>app</IT></SUB> (8)
where C is the same as in model 1. The ionic currents INa, IK, IL, Itonic-e, and Iapp are identical to those in model 1. INaP is identical to INaP-h in model 1 except that it does not possess the inactivation term h.
<IT>I</IT><SUB><IT>KS</IT></SUB><IT>=</IT><IT><A><AC>g</AC><AC>&cjs1171;</AC></A></IT><SUB><IT>KS</IT></SUB><IT>k</IT>(<IT>V</IT><IT>−</IT><IT>E</IT><SUB><IT>K</IT></SUB>) (9)
where <A><AC>g</AC><AC>&cjs1171;</AC></A>KS = 5.6 nS, theta k = -38 mV, sigma k = -6 mV, <A><AC>&tgr;</AC><AC>&cjs1171;</AC></A>k = 10,000 ms, and EK = -85 mV. The voltage dependency of the gating kinetics for k was chosen so that IKS was activated by the firing of action potentials and was also dynamically active in the subthreshold potential range from -60 to -45 mV. Figure 2C1 illustrates kinfinity (V) and tau k(V). Figure 2C2 illustrates INaP + IKS when k = kinfinity (V) (slow voltage ramp with other membrane currents eliminated) and k = 0 (removal of IKS, similar to a fast voltage ramp). The shapes of these plots of INaP + IKS are similar to those of INaP-h illustrated in Fig. 2B2.

The slow variables h in model 1 and k in 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 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.


    RESULTS
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

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 changing Iapp to control the neuron's level of depolarization is equivalent to varying EL (Iapp = <A><AC>g</AC><AC>&cjs1171;</AC></A>LDelta EL).



View larger version (37K):
[in this window]
[in a new window]
 
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 (Iapp), changing [K+]o (which in turn affects EK) or adjusting the level of tonic synaptic input, the mean level of which is represented by gtonic-e.

Figure 4 illustrates the effect of varying EL on the oscillatory activity of model 1. The time course of membrane potential is illustrated in Fig. 4A, 1-4, and the time course of the INaP-h inactivation (h) and instantaneous membrane conductance are illustrated in Fig. 4B, 1-4. As EL 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 as EL is increased. Vmin, the minimum membrane potential encountered during the silent phase of the burst cycle, varies from -58 mV at the onset of bursting activity (EL = -60.5 mV) to -48 mV (EL = -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, INaP-h is inactivated progressively with the firing of each action potential (Fig. 4D, 2 and 3). During the interburst phase, inactivation is removed and the instantaneous membrane conductance gradually increases as INaP-h recovers from inactivation. Similar results to those in Fig. 4A are obtained if EL is separated into K+ and Na+ components and EK is varied.



View larger version (45K):
[in this window]
[in a new window]
 
Fig. 4. Dynamic response of model 1 as a function of EL. Left: (A, 1-4) membrane potential. Middle: (B, 1-4) instantaneous membrane conductance (gm) and inactivation of INaP-h (h) as cell is depolarized (EL = -65, -60, -57.5, and -54 mV for 1-4, respectively). C2/D2 and C3/D3: expanded time course of V and h and the spike-frequency profile (open circle ) 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 varying gtonic-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 of gtonic-e on the mode of activity and pattern of bursting are similar to those of EL (quantified in the following text).

Figure 5 illustrates the effect of varying EL on the oscillatory activity of model 2. The time course of membrane potential is illustrated in Fig. 5A, 1-4, and the time course of the IKS 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 EL is increased, the model progresses through regimes of silence, bursting, and beating activity. During oscillatory bursting, the burst is initiated by INaP and terminates when IKS has been activated sufficiently. The spike frequency decays throughout the active phase of the burst. The current IKS is activated slowly and progressively with the firing of each action potential and deactivates during the interburst interval. The instantaneous membrane conductance (gm) decreases throughout most of the silent phase and increases toward the end of the phase. This differs from model 1, where gm steadily increases throughout the silent phase. This difference in the time course of gm may provide a means for differentiating between the mechanisms of the two models in vitro (see DISCUSSION).



View larger version (43K):
[in this window]
[in a new window]
 
Fig. 5. Dynamic response of model 2 as a function of EL. Left: (A, 1-4) membrane potential. Right: (B, 1-4) instantaneous membrane conductance (g) and activation of IKS (k) as cell is depolarized (EL = -65, -59.5, -50, and -40 mV for 1-4, respectively). C2/D2 and C3/D3: expanded time course of V and k and the spike-frequency profile (open circle ) 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 IKS and INaP. 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 gm did 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 with EL or gtonic-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 <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP. The burst duration, burst period, and frequency of bursting and beating are plotted as EL (Fig. 6A) and gtonic-e (Fig. 6B) are varied. Equivalent values of Iapp with EL fixed to -65 mV are also shown on A. 



View larger version (26K):
[in this window]
[in a new window]
 
Fig. 6. Intrinsic and extrinsic mechanisms for frequency control in model 1. Panels show effect of varying EL (A) and gtonic-e (B) on burst duration, burst period, burst frequency, and beat frequency. These curves were calculated for multiple values of <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP. Bursting activity does not occur for <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP = 2.0 nS in A and for <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP = 2.0 and 2.4 nS in B. IAPP below EL scale shows equivalent variation in IAPP with EL fixed to -65 mV.

In both panels, there is a minimum value of <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP for which oscillatory bursting can occur. If <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP is too low, the only supported modes of activity are quiescence and beating as shown for model 1 in Fig. 7. For all values of gNaP 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 <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP increases the dynamic range of values of EL where bursting is supported. We also have observed (not shown) that increasing <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP increases the range of values of Vmin encountered as a function of EL during oscillatory bursting. At a given level of EL, increasing <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP increases burst frequency.



View larger version (41K):
[in this window]
[in a new window]
 
Fig. 7. Intrinsic modes of single-cell activity of model 1 as a function of <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP and EL. When <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP exceeds 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 EL. Range of <A><AC>g</AC><AC>&cjs1171;</AC></A>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 of EL on model 2 bursting for a range of values of the subthreshold conductances <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP and <A><AC>g</AC><AC>&cjs1171;</AC></A>KS. For all values of <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP and <A><AC>g</AC><AC>&cjs1171;</AC></A>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 by EL 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 of EL (and Iapp) approximately twice as large as the range of EL where bursting occurred in model 1. These trends persisted as the other parameters of the subthreshold currents (e.g., theta k and sigma k) were varied. Similar results were obtained when gtonic-e was varied.



View larger version (28K):
[in this window]
[in a new window]
 
Fig. 8. Mechanisms for frequency control in model 2. Panels show effect of varying EL on burst duration, burst period, burst frequency, and beat frequency of model 2. These curves were calculated for multiple values of <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP (A) and <A><AC>g</AC><AC>&cjs1171;</AC></A>KS (B). IAPP below EL scale shows equivalent variation in IAPP with EL fixed to -65 mV.

The large range of Iapp in model 2 may be attributed to the difference in whole cell conductance (gm) between the two models. Model 2 has an additional conductance, <A><AC>g</AC><AC>&cjs1171;</AC></A>KS, which is not present in model 1. This conductance is activated during bursting and beating. The increased gm means that the "gain" of the cell is decreased. A larger increment in Iapp 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 larger Iapp 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 <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP and <A><AC>g</AC><AC>&cjs1171;</AC></A>KS. In both models (1 and 2), the critical value of EL or Iapp that changes the cell from quiescence to bursting is quite sensitive to <A><AC>g</AC><AC>&cjs1171;</AC></A>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 <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP. In model 2, <A><AC>g</AC><AC>&cjs1171;</AC></A>KS affects this critical level only modestly, partly because the driving force for <A><AC>g</AC><AC>&cjs1171;</AC></A>KS is much lower than that of <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP and partly because it possesses slow (compared with INaP) 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 <A><AC>g</AC><AC>&cjs1171;</AC></A>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 EL and Iapp are almost independent of <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP. The transition is affected by <A><AC>g</AC><AC>&cjs1171;</AC></A>KS because this current controls burst termination (Fig. 5B3). Although these sensitivities for the silence to bursting transition can be understood qualitatively by considering the subthreshold currents (Isub, see following text) alone, the corresponding sensitivities for the bursting to beating transition are difficult to intuit because the components of Isub 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 varying EL, 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 (Isub triple-bond  IL + INaP-h for model 1, Isub triple-bond  IL +INaP + IKS for model 2). This explanation suffices because without the action potential currents INa and IK, reduced models having only the Isub currents display oscillatory and nonoscillatory states (hyperpolarized silence, oscillatory activity, depolarized silence) that are qualitatively similar to those in the corresponding full models. Figure 9 illustrates I-V curves corresponding to Isub in cases shown in Figs. 4 and 5. The SS curves (thick lines) represent Isub versus V when the slow negative feedback process (h inactivation for model 1, k activation for model 2) is set to steady state [(i.e., h = hinfinity (V) in model 1, and k = kinfinity (V) in model 2]). The QSS curves (thin lines) represent Isub 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. 9A, 1-4). An analogous explanation suffices for model 2.



View larger version (28K):
[in this window]
[in a new window]
 
Fig. 9. Subthreshold currents of model 1 (INaP-h + IL) (left) and model 2 (INaP + IKS + IL) (right) as EL is varied. Values of EL correspond to values used in Figs. 4 and 5 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 = hinfinity (V) in A, 1-4, k = kinfinity (V) in B, 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. In A1 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 and II, top). In A4 and B4, the QSS I-V curve was generated by setting h or k to the mean value during tonic firing (see Figs. 4B4 and 5B4).

For all values of EL the SS I-V curves 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 (CM/gm), 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-V curves provide any information on the dynamics of h itself.

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

For EL equal to -60 mV (Fig. 4A2) or -57.5 mV (Fig. 4A3), 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), Isub(t), and QSS I-V(V, t), is available on the World-Wide Web at http://intra.ninds.nih.gov/smith/movie/moviejnp-99.html. Figure 9A, 2 and 3, shows Isub at the maximum (a) and minimum (b) values of h during one complete burst cycle. The burst begins at a value of h (a) where Isub 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 INaP-h to Isub, shifting the QSS I-V curve outward. Firing persists until Isub is net outward at Vmin, 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-V curve (b). As h gradually deinactivates, the QSS I-V curve moves inward, depolarizing the cell as the pseudosteady equilibrium point drifts rightward. Another burst begins when Isub 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 to A3, Fig. 4), the values of h encountered during one burst cycle are biased toward smaller values, i.e., less inward current from INaP-h is necessary to counterbalance the reduced outward current from shifting EL. In addition, the dynamic range of h is reduced, from a Delta h of ~0.1 in Fig. 4A2 to <0.02 in Fig. 4A3. This is the essence of the voltage-dependent frequency control in the model: as the cell is depolarized toward higher values of EL, the cell does not hyperpolarize as far at the end of each burst cycle, and less time is necessary to recover from inactivation of INaP-h to begin the following burst. At all values of EL where bursting occurs, INaP-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. 4A4, EL = -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. 9A4, generated for this mean value of h (0.315), illustrates that on average Isub 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 <A><AC>g</AC><AC>&cjs1171;</AC></A>Na = 0. This predicts what we expect in experiments when INa is blocked (e.g., with QX-314 or TTX-see DISCUSSION), 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. 4A, 1-4, and 5A, 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.



View larger version (17K):
[in this window]
[in a new window]
 
Fig. 10. Oscillatory model behavior in the absence of fast Na+ currents. Panels show effects of setting <A><AC>g</AC><AC>&cjs1171;</AC></A>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. 4A, 1-4, and 5 A, 1-4, respectively. Parameters for A5 and B5 are EL = -53 mV and EL = -32 mV, respectively. - - -, -60 and -30 mV references.

Cells in hyperpolarized silent mode may retain burst excitability if <A><AC>g</AC><AC>&cjs1171;</AC></A>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. 11A). Similar responses occur across a wide range of values of EL and <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP where the model is initially silent, providing that <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP is not too small. At rest, INaP-h is relatively deactivated (m is low) and deinactivated (h is high). A brief depolarization activates INaP-h. Because INaP-h inactivates slowly, INaP-h remains relatively deinactivated and a larger conductance for INaP-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.



View larger version (19K):
[in this window]
[in a new window]
 
Fig. 11. Transient burst responses. A: brief but sufficiently large transient depolarizing input can elicit a single sustained burst. Model was at rest (EL = -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 (EL = -62 mV), and a hyperpolarizing current pulse (500 ms, 60 pA) was applied. In both cases, the parameters EL and <A><AC>g</AC><AC>&cjs1171;</AC></A>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. 11B). A sustained hyperpolarizing input removes some inactivation from INaP-h (i.e., h increases). On abrupt release of the hyperpolarizing input, the model depolarizes. Because inactivation of INaP-h develops slowly, INaP-h has an increased conductance as V 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 of INaP-h is a slow process), the hyperpolarizing stimulus is of sufficient magnitude (so that h increases sufficiently), and the resting neuron is not too hyperpolarized (so that sufficient h remains to be deinactivated). For example, the response of Fig. 11B was generated for EL = -62 mV. When EL = -65 mV, the resting value of h is sufficiently large (0.92) that even a large hyperpolarization (to bring h nearly to 1) does not produce an increase in the conductance of INaP-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 (b1) terminates the burst, and the following interburst interval is significantly shortened, with the next burst firing at b2. A transient current pulse applied late in the burst (c1) also terminates the burst, with the next burst firing at c2. 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).



View larger version (22K):
[in this window]
[in a new window]
 
Fig. 12. Resetting response. Brief but sufficiently strong transient hyperpolarizing input can reset a burst. Simulations performed with EL = -59 mV, which corresponds to a burst period of ~4 s. Control period shown in A, with bursts occurring at times alpha 1 and alpha 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.


    DISCUSSION
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

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 (<A><AC>g</AC><AC>&cjs1171;</AC></A>NaP) and baseline membrane potentials (EL); 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.

BURST INITIATION. 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 in Helix (Swandulla and Lux 1985) and Aplysia (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 IH 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 IH 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 IH 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 IH 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 (IR). IR 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 (INaP). INaP 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 of INaP has been quantitatively characterized in mouse and rat neocortical neurons (Brown et al. 1994; Fleidervish et al. 1996; French et al. 1990). We adopted INaP as the current responsible for burst initiation in our model. Ongoing experiments in our laboratory seek to identify the existence of INaP 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 INaP: theta m = -50 and -48 mV, respectively, and absolute values of sigma m in the range of 4-9 mV. Our choice of theta m (-40 mV) is depolarized from these values. However, similar voltage-dependent frequency control (at more hyperpolarized potentials) is obtained from the model if theta m is reduced to -48 mV with a similar shift in theta h. The formulation of minfinity (V) for INaP in our model is consistent with observations from cortical neurons where the threshold for INaP activation is >= 10 mV hyperpolarized to the threshold for INa 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 INaP dominates the excitability of the membrane at subthreshold potentials (Stafstrom et al. 1982). This is evident when comparing m(V) for INa with minfinity (V) for INaP in Fig. 2 (A and B1).

BURST TERMINATION. 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 beta -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 of INaP-h have not been as well characterized; the parameters used in our models are consistent with available measurements. Our hinfinity (V) is similar in shape to the voltage dependency of the slow inactivation of INaP-h reported by Fleidervish et al. (1996). The steady-state I-V plots of INaP-h in our model (Fig. 2B2) 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 INaP 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 INaP-h as 2.06 s at +20 mV and the time constant for recovery from slow inactivation of INaP-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 tau h(V) (Fig. 2B1), although they lie at extreme ends of the bell-shaped curve tau h(V). We used a value of 10 s for <A><AC>&tgr;</AC><AC>&cjs1171;</AC></A>h to produce burst periods on the time scales that are observed experimentally.

The hypothetical burst termination mechanism in model 2 requires a slowly activating voltage-dependent K+ current. Slow voltage-gated Phe-Met-Arg-Phe-NH2 (FMRF-amide) sensitive K+-currents (IKF) 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 Iapp or gtonic-e, or varying EL (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: <A><AC>g</AC><AC>&cjs1171;</AC></A>KS, <A><AC>g</AC><AC>&cjs1171;</AC></A>Nap, theta m, sigma m, theta k, sigma k) in model 2 and have not identified parameter regimes where model 2 exhibits such behavior.

We found that by varying additional parameters of IKS we can mathematically transform the ionic current equations of model 1 into a "nearly dynamically equivalent" model that possess the ionic current equations of model 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 IKS -55 mV or incorporate an additional fast-activating gating variable to IKS such that IKS is only active at membrane potentials above -55 mV. With either of these manipulations to model 2, gm decreases 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: Isub approx  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: Delta V/Delta h or Delta V/Delta k. For model 1, we find that Delta V/Delta h is proportional to INaP (i.e., INaP-h with h = 1), while for model 2 Delta V/Delta k is proportional to <A><AC>g</AC><AC>&cjs1171;</AC></A>KS(V - EK) (i.e., IKS with k = 1). Because INaP-h is almost totally deactivated (m approx  0) in the voltage range of the interburst phase we see that Delta V/Delta h is much less than Delta V/Delta k. 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 Iapp in model 1 than in model 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 by model 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 (see Rekling 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 present model 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 INaP-h recovers from inactivation (Fig. 4B, 1-4). The recovery of h and the gradual depolarization of the membrane potential (via the fast activation of INaP-h) both contribute to an increase in <A><AC>g</AC><AC>&cjs1171;</AC></A>Nap-h, and thus gm. The presence of additional subthreshold currents with voltage-dependent conductances (e.g., IR or IH), however, could alter this trend.

In model 2 during the silent phase of the burst cycle, there are two competing effects with respect to gm. The deactivation of IKS contributes toward a reduction of gm, while the membrane's gradual depolarization activates INaP, contributing toward an increase in gm. Thus the overall change in gm in model 2 depends on the relative strength of INaP and IKS. In Fig. 5 it is shown that gm 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 IKS, and on hyperpolarization (as long as we hyperpolarize above EK) 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 INaP-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 6B 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 EL or Iapp. Similar results to Fig. 6 are obtained as EK is varied (with IL 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 as EK or gtonic-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. INaP, 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 INaP 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. IH and IT 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 that INaP-h is another current mechanism promoting posthyperpolarization bursting (Fig. 11B). 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 Iapp is increased. Similar effects of Iapp 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 beta -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-beta 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 Iapp. Thus we reexamined a few of these beta -cell models (Sherman and Rinzel 1992; Sherman et al. 1990; Smolen and Keizer 1992) to study the effects of Iapp 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 Iapp, 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 in Smith 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 <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP may be viewed as a tunable parameter that scales the dynamic range of bursting (Figs. 6 and 7). As <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP is increased, two effects are observed: the range of values of EL where bursting is supported is increased and the range of burst periods possible as EL is varied is increased. The burst duration is primarily a function of EL and relatively independent of values of <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP.

Model 2 offers two distinct subthreshold conductances for controlling the model's burst dynamics, <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP and <A><AC>g</AC><AC>&cjs1171;</AC></A>KS (Fig. 8). Decreasing <A><AC>g</AC><AC>&cjs1171;</AC></A>KS and/or increasing <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP increases burst duration. <A><AC>g</AC><AC>&cjs1171;</AC></A>KS and <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP have complementary roles with respect to controlling the frequency of bursting: <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP controls the minimal value of EL where bursting is supported and the maximum burst period, whereas <A><AC>g</AC><AC>&cjs1171;</AC></A>KS controls the maximal value of EL (or gtonic-e or Iapp) where bursting is supported and the minimal burst period. Although <A><AC>g</AC><AC>&cjs1171;</AC></A>KS appears to set an upper bound on the range of EL values that support bursting, Figure 8B2 illustrates that <A><AC>g</AC><AC>&cjs1171;</AC></A>KS has a minimal effect on the overall burst period, which may be considered to be primarily a function of EL. Thus <A><AC>g</AC><AC>&cjs1171;</AC></A>KS may serve as a mechanism for varying duty cycle with minimal effects on the overall burst frequency.

In both models, gL 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 gL results in quiescence. As gL is reduced, bursting occurs at a frequency dependent on the value of gL. Sufficiently low values of gL 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 of INaP: 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 <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP determines whether the neuron will exhibit oscillatory bursting: at low values of <A><AC>g</AC><AC>&cjs1171;</AC></A>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 <A><AC>g</AC><AC>&cjs1171;</AC></A>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 <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP and EL 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 of INaP and synaptic interactions, synchronous oscillations can occur in the coupled population even if none (low INaP) or few of the cells exhibit oscillatory bursting. As discussed earlier, <A><AC>g</AC><AC>&cjs1171;</AC></A>NaP and gL 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.

Summary

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.


    ACKNOWLEDGMENTS

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.


    FOOTNOTES

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.

Received 23 September 1998; accepted in final form 9 February 1999.


    REFERENCES
TOP
ABSTRACT
INTRODUCTION
METHODS
RESULTS
DISCUSSION
REFERENCES

0022-3077/99 $5.00 Copyright © 1999 The American Physiological Society