Dual oscillator model of the respiratory neuronal network generating quantal slowing of respiratory rhythm

Abstract

We developed a dual oscillator model to facilitate the understanding of dynamic interactions between the parafacial respiratory group (pFRG) and the preBötzinger complex (preBötC) neurons in the respiratory rhythm generation. Both neuronal groups were modeled as groups of 81 interconnected pacemaker neurons; the bursting cell model described by Butera and others [model 1 in Butera et al. (J Neurophysiol 81:382–397, 1999a)] were used to model the pacemaker neurons. We assumed (1) both pFRG and preBötC networks are rhythm generators, (2) preBötC receives excitatory inputs from pFRG, and pFRG receives inhibitory inputs from preBötC, and (3) persistent Na⁺ current conductance and synaptic current conductances are randomly distributed within each population. Our model could reproduce 1:1 coupling of bursting rhythms between pFRG and preBötC with the characteristic biphasic firing pattern of pFRG neurons, i.e., firings during pre-inspiratory and post-inspiratory phases. Compatible with experimental results, the model predicted the changes in firing pattern of pFRG neurons from biphasic expiratory to monophasic inspiratory, synchronous with preBötC neurons. Quantal slowing, a phenomena of prolonged respiratory period that jumps non-deterministically to integer multiples of the control period, was observed when the excitability of preBötC network decreased while strengths of synaptic connections between the two groups remained unchanged, suggesting that, in contrast to the earlier suggestions (Mellen et al., Neuron 37:821–826, 2003; Wittmeier et al., Proc Natl Acad Sci USA 105(46):18000–18005, 2008), quantal slowing could occur without suppressed or stochastic excitatory synaptic transmission. With a reduced excitability of preBötC network, the breakdown of synchronous bursting of preBötC neurons was predicted by simulation. We suggest that quantal slowing could result from a breakdown of synchronized bursting within the preBötC.

1 Introduction

Two putative respiratory rhythm generators, the parafacial respiratory group (pFRG) (Onimaru et al. 1988; Onimaru and Homma 2003) and the preBötzinger complex (preBötC) (Smith et al. 1991), exist in the mammalian brainstem. The two neuronal groups, separately identified, are thought to play important roles in the generation and maintenance of respiratory rhythm (Feldman and Janczewski 2006; Onimaru and Homma 2006; Oku et al. 2007) at least in the neonatal period, although the primary site of the respiratory rhythm generation is still controversial (Feldman and Del Negro 2006). There are a number of neuronal behaviors that characterize the respiratory network in brainstem spinal cord preparations of neonatal rodents. First, the firing pattern of preinspiratory (Pre-I) neurons, a subset of neurons within pFRG, typically consists of three phases: (1) firing preceding the inspiratory phase, (2) suppression of firing during the inspiratory phase, originally referred to as inspiration-related inhibition of the Pre-I firing (IIPI), and (3) firing during the post-inspiratory phase, which is not seen in preBötC preinspiratory neurons (Smith et al. 2007). Second, IIPI consistently disappears after removing chloride-mediated inhibition, resulting in overlapping of Pre-I neuronal firing and inspiratory neuronal firing (Onimaru et al. 1990). Third, opioids do not affect Pre-I neurons but depress preBötC inspiratory neurons (Gray et al. 1999). As a consequence of dynamic interactions between pFRG and preBötC, respiratory periods jump non-deterministically to integer multiples of the control period, and the phenomenon was referred to as quantal slowing (Mellen et al. 2003; Barnes et al. 2007). Synaptic connections and interactions between the neuronal groups to explain the characteristic behaviors have been proposed (Onimaru et al. 1990; Ballanyi et al. 1999). Specifically, the observation of quantal slowing has led to a hypothesis that two rhythmically active networks interact to generate respiratory rhythm and quantal slowing results from transmission failure from pFRG to preBötC networks (Mellen et al. 2003). However, in such a complex dynamic system, it is imperative to evaluate descriptive assumptions by simulating phenomena using a computational model.

In the present study, we demonstrate a computational model consisting of two groups of bursting neurons, representing pFRG and preBötC neuronal groups. We assume that pFRG provides excitatory inputs to preBötC and preBötC provides inhibitory connections to pFRG (Onimaru et al. 1990; Ballanyi et al. 1999). We then analyzed the influences of changes in strengths of excitatory and inhibitory synaptic connections between the two neuronal groups on the coupling pattern of bursting activities in these groups. Finally, we evaluated the effects of reducing excitability of preBötC neurons on the respiratory period, and considered the mechanism that causes non-deterministic jump of respiratory period to integer multiples of the control periods. Interestingly, we found that the quantal slowing occurs without assuming suppressed excitatory synaptic transmission.

Similar attempts to simulate quantal slowing phenomenon by dynamic interactions between the two rhythm generating networks has been made by Joseph and Butera (2005) and Wittmeier et al. (2008) earlier. In both studies, pFRG and preBötC were represented by only a pair of pacemaker neurons, presuming that any ‘lumped’ effects of the interaction should apply at the full-blown network level (Wittmeier et al. 2008). Joseph and Butera (2005) used a simple canonical model of pacemaker neurons to represent the pFRG and preBötC. However, it is difficult to correlate cellular properties with model parameters in such a canonical model. Further, in their result, quantal slowing accompanied abrupt changes in firing pattern of Pre-I neurons, i.e., emergence of post-inspiratory burst that was absent in the control period, which is inconsistent with experimental findings (Mellen et al. 2003). Wittmeier et al. (2008) adopted a Hodgkin-Huxley type realistic neuronal model described by Butera et al. (1999a) and substituted stochastic synaptic transmission for the deterministic excitatory synaptic inputs to preBötC from pFRG to simulate quantal slowing. However, there is no physiological basis supporting that opioids cause stochastic synaptic transmission. Moreover, if the stochastic nature of synaptic transmission is essential for producing quantal slowing then the phenomenon is not entirely a characteristic of the model under consideration but also of a undefined additional mechanism resulting in stochastic dynamics of the synaptic transmission. In the present study, we have conducted a large scale simulation by modeling both pFRG and preBötC as groups of 81 pacemaker neurons to test whether the essential features of the activity of a rhythm generating network may be adequately captured by replacing it by a single pacemaker neuron, and more specifically, to test whether the quantal slowing phenomenon requires stochastic dynamics of the synaptic transmission. Our simulation predicted that the reduction of excitability in preBötC network could result in intermittent failure of preBötC neurons to produce synchronized bursting and could cause quantal slowing even if we did not assume stochastic dynamics of the synaptic transmission. We suggest that the quantal slowing phenomenon is an intrinsic characteristic of the rhythm-generating network and stochastic synaptic transmission is unnecessary to explain the quantal slowing phenomenon. Our results are in agreement with the general conception that quantal slowing phenomenon results due to general loss of sensitivity within the preBötC group but, more importantly, they suggest towards the possibility of a different mechanism for its cause.

2 Methods

2.1 Formulation of dual oscillator model

The model consists of two neuronal groups of bursting neurons, NeuronGroup1 and NeuronGroup2; representing the pFRG and preBötC rhythm generating networks, respectively. We assume that both pFRG and preBötC independently constitute rhythm-generating networks, and pFRG network excites preBötC network, and preBötC network inhibits pFRG network (Fig. 1(a)). In our model, we assume that bursting neurons within each of NeuronGroup1 and NeuronGroup2 are mutually coupled by excitatory synaptic connections to form a rhythm generating network. Therefore, the constituent neurons of NeuronGroup2 are excitatory neurons. It is unrealistic to assume that the excitatory neurons of NeuronGroup2 inhibit NeuronGroup1. A more realistic model would be a three neuronal group model as depicted in Fig. 1(b); the third neuronal group, NeuronGroup3, is formed of inhibitory neurons, which receives excitatory synaptic input from NeuronGroup2 and provides inhibitory synaptic input to NeuronGroup1 (Kuwana et al. 2006). However, for simplicity, if we assume that the bursting dynamics of NeuronGroup3 follow the bursting activities of NeuronGroup2, then the NeuronGroup3 in Fig. 1(b) essentially becomes a copy of NeuronGroup2. Thus, the need of explicitly incorporating NeuronGroup3 is eliminated and computational effort required is minimized. Therefore, we use the model depicted in Fig. 1(a) as our Dual Oscillator Model.

Fig. 1

(a) Dual Oscillator Model: NeuronGroup1 (pFRG) provides excitatory synaptic input to NeuronGroup2 (preBötC) and NeuronGroup2 provides inhibitory synaptic input directly to NeuronGroup1. (b) Three neuronal group model: Similar to Dual Oscillator Model except that an inspiratory interneuron group, NeuronGroup3 is added. NeuronGroup3 receives excitatory inputs from NeuronGroup2 and inhibit NeuronGroup1. Open and closed circles represent excitatory and inhibitory synaptic connections, respectively

2.2 Formulation of neuronal groups

Both neuronal groups, NeuronGroup1 and NeuronGroup2, were modeled as groups of 81 interconnected pacemaker neurons. The mathematical model of bursting cells described by Butera and others (model 1 in Butera et al. 1999a) were used to model the pacemaker neurons in our model. We assume that each cell receives excitatory and inhibitory synaptic inputs. Therefore, the rate of change in membrane potential (V) of a cell is modified from the original model as:

$$ C\frac{dV}{dt}=-I_{\rm NaP}-I_{\rm Na}-I_{\rm K}-I_{\rm L}-I_{\rm syn(e)}-I_{\rm syn(i)} $$

(1)

where C is the whole cell capacitance (pF), V is membrane potential (mV), t is time (ms), I _NaP is a persistent Na⁺ current, I _Na is a fast Na⁺ current, I _K is a delayed-rectifier K⁺ current, I _L is a passive leakage current, I _syn(e) is the sum of excitatory synaptic currents, and I _syn(i) is the sum of inhibitory synaptic currents. C is set to 21 pF (Smith et al. 1992; Butera et al. 1999a). Equations and parameter sets to compute intrinsic membrane currents (I _NaP, I _Na, I _K and I _L) are those presented in the previous paper (Butera et al. 1999a) unless otherwise noted.

We assume that each neuron receives excitatory synaptic inputs from other neurons belonging to the same group L as well as from neurons belonging to a different group M:

$$\begin{array}{rll} I_{\rm syn(e)}&=&\left(\sum\limits_{l\in L}g^{l}_{\rm syn(int)}\cdot s^{l}_{\rm int}\right)\left(V-E_{\rm syn(e)}\right) \\&&+\left(\sum\limits_{m\in M}g^{m}_{\rm syn(ext)}\cdot s^{m}_{\rm ext}\right)\left(V-E_{\rm syn(e)}\right) \end{array}$$

(2)

where E _syn(e) (= 0 mV) is the reversal potential of non-NMDA excitatory amino-acid receptors. $g^{l}_{\rm syn(int)}$ and $g^{m}_{\rm syn(ext)}$ are the synaptic conductances for input from neurons l ∈ L and m ∈ M, respectively. $s^{l}_{\rm int}$ and $s^{m}_{\rm ext}$ are synaptic gating variables for these synaptic inputs. To reduce computational costs, we assumed that synaptic conductances $g^{l}_{\rm syn(int)}$ and $g^{m}_{\rm syn(ext)}$ in Eq. (2) are uncorrelated to the corresponding synaptic gating variables $s^{l}_{\rm int}$ and $s^{m}_{\rm ext}$, respectively. Consequently, we simplify Eq. (2) as:

$$ I_{\rm syn(e)}=\left(\bar{g}_{\rm syn(int)}\cdot\bar{s}_{\rm int}+\bar{g}_{\rm syn(ext)}\cdot\bar{s}_{\rm ext}\right)\left(V-E_{\rm syn(e)}\right) $$

(3)

where $\bar{g}_{\rm syn(int)}$ and $\bar{g}_{\rm syn(ext)}$ are overall synaptic conductance for inputs from neurons l ∈ L and m ∈ M, respectively. The values of $\bar{g}_{\rm syn(int)}$ and $\bar{g}_{\rm syn(ext)}$ were randomly assigned from a uniformly distributed probability density function ranging between 0 and 10 nS and between 0 and $\bar{g}_{\rm syn(ext)\underline{~}max}$, respectively. The rationale for using uniformly distributed probability density function for synaptic conductance is discussed in Section 4.3. $\bar{s}_{\rm int}$ and $\bar{s}_{\rm ext}$ are mean synaptic gating variables for these synaptic inputs, which are expressed as:

$$ \bar{s}_{\rm int}=\overline{\sum\limits_{l\in L}s^{l}_{\rm int}}~,\bar{s}_{\rm ext}=\overline{\sum\limits_{m\in M}s^{m}_{\rm ext}} $$

(4)

Similarly, synaptic inputs from an inhibitory neuronal group N are expressed as:

$$ I_{\rm syn(i)}=\bar{g}_{\rm syn(inh)}\cdot\bar{s}_{\rm inh}\left(V-E_{\rm Cl}\right) $$

(5)

where $\bar{g}_{\rm syn(inh)}$ is the overall synaptic conductance for inputs from neurons n ∈ N, $\bar{s}_{\rm inh}$ are the mean synaptic gating variable for inhibitory inputs, and E _Cl (= − 90 mV) is the equilibrium potential of Cl⁻. The values of $\bar{g}_{\rm syn(inh)}$ were randomly assigned from a uniformly distributed probability density function ranging between 0 and $\bar{g}_{\rm syn(inh)\underline{~}max}$.

Presynaptic action potentials activate synaptic gating variables. We applied the same, fast (decay time constant τ _s  = 5 ms) kinetics to all synaptic gating variables (Butera et al. 1999b):

$$ \frac{ds_{i}}{dt}=\frac{\left[\left(1-s_{i}\right)\cdot s_{\infty}\left(V_{i}\right)-s_{i}\right]}{\tau_{s}} $$

(6)

$s_{\infty}\left(V_{i}\right)$ determines the steady-state postsynaptic receptor activation based on the membrane potential of presynaptic neuron i:

$$ s_{\infty}\left(V_{i}\right)=\frac{1}{1+\exp\left[\frac{\left(V_{i}-\theta\right)}{\sigma}\right]} $$

(7)

where θ = − 10 mV and σ = − 5 mV.

In the present model, I _NaP produces a slowly depolarizing current and endows cells with a bursting property. The inherent burst frequency of a single bursting neuron is changed by the persistent Na⁺ current conductance, $\bar{g}_{\rm NaP}$, with a fixed value of the reversal potential of passive leakage current E _L (− 59 mV). Bursting activity does not occur for $\bar{g}_{\rm NaP}\leq$ ≤ 2.45 nS (Butera et al. 1999a). In order to simulate an inhomogeneous distribution of inherent burst frequency within a group, the values of $\bar{g}_{\rm NaP}$ were randomly assigned from a uniformly distributed probability density function ranging between ($\bar{g}_{\rm NaP\underline{~}max}-0.5$ nS) and $\bar{g}_{\rm NaP\underline{~}max}$. Consequently, the inherent burst frequency of the neuronal group depends on $\bar{g}_{\rm NaP\underline{~}max}$. The firing of Pre-I neurons always precedes the inspiratory burst recorded from the fourth cervical ventral root of the spinal cord (C4VR), although Pre-I burst sometimes does not accompany the inspiratory C4VR burst (Okada et al. 2007). In order to simulate this situation, the inherent burst frequency of NeuronGroup1 must be greater than that of NeuronGroup2. Subsequently, in the present model, we set $\bar{g}_{\rm NaP\underline{~}max}=4.0$ nS for NeuronGroup1 and $\bar{g}_{\rm NaP\underline{~}max}=3.0$ nS for NeuronGroup2. For these values of $\bar{g}_{\rm NaP\underline{~}max}$ for NeuronGroup1 and NeuronGroup2, their interburst period was, respectively, found to be 2.69 ± 0.13 s and 4.99 ± 0.42 s (computed from 9 simulation results).

We set $\bar{g}_{\rm syn(ext)} = 0$ nS for NeuronGroup1 and $\bar{g}_{\rm syn(inh)} = 0$ nS for NeuronGroup2. $\bar{s}_{\rm inh}$ of NeuronGroup1 was calculated from membrane potentials of NeuronGroup2 neurons and $\bar{s}_{\rm ext}$ of NeuronGroup2 was calculated from membrane potentials of NeuronGroup1 neurons. To evaluate the effects of changes in strengths of excitatory and inhibitory synaptic connections on the coupling pattern between pFRG and preBötC bursting activities, we varied $\bar{g}_{\rm syn(inh)\underline{~}max}$ of NeuronGroup1 ranging between 0 and 6 nS and $\bar{g}_{\rm syn(ext)\underline{~}max}$ of NeuronGroup2 ranging between 0 and 2 nS.

2.3 Computational details

The differential equations were solved numerically using the fourth-order Runge-Kutta equation with a step size of 0.05 ms. However, the neuronal states were saved at an interval of 0.5 ms, which was sufficient to capture every action potential spikes of bursting neurons, to keep the stored files to manageable size. All algorithms were implemented in double precision routine in C$\sharp$.NET language and run on Pentium-based Windows XP computers. We waited for initial 60 s of simulation time for the behavior of the network to stabilize and then the simulated data were saved.

3 Results

As mentioned earlier, in the dual oscillator model, NeuronGroup1 and NeuronGroup2 represents pFRG and preBötC, respectively. However, during the presentation of results and discussion, we restrict the terms NeuronGroup1 and NeuronGroup2 to refer to the modeled pFRG and preBötC, respectively. Further, throughout the present model study, we define the period where preBötC/NeuronGroup2 neurons burst synchronously as the inspiratory phase, and the rest of period where a majority of preBötC/NeuronGroup2 neurons are quiescent as the expiratory phase.

3.1 Coupling pattern between pFRG and preBötC

We observed different patterns of coupling between NeuronGroup1 and NeuronGroup2 depending on strengths of excitatory and inhibitory synaptic connections ($\bar{g}_{\rm syn(ext)\underline{~}max}$ and $\bar{g}_{\rm syn(inh)\underline{~}max}$, respectively). Figure 2(a) depicts typical raster plots for various coupling modes observed. Horizontal segments in Fig. 2(a) represent bursting periods of all the 2 × 81 neurons constituting the NeuronGroup1 and NeuronGroup2. For each burst of each neuron, the bursting period is defined as the period beginning from the instant when the membrane potential of the neuron rises above − 20 mV while delivering the first spike of the burst to the instant when the membrane potential falls below − 20 mV while delivering the final spike of the burst. Membrane potential trajectories of a typical neuron within NeuronGroup1 and NeuronGroup2 for the various coupling modes are depicted in Fig. 2(c); averaged population activities of NeuronGroup1 and NeuronGroup2 are depicted in Fig. 2(b).

Fig. 2

(a) Raster plots of various coupling modes—biphasic ($\bar{g}_{\rm syn(ext)\underline{~}max}=0.8$ nS, $\bar{g}_{\rm syn(inh)\underline{~}max}=4.5$ nS), monophasic ($\bar{g}_{\rm syn(ext)\underline{~}max}=2.0$ nS, $\bar{g}_{\rm syn(inh)\underline{~}max}=5.0$ nS), Synchronous ($\bar{g}_{\rm syn(ext)\underline{~}max}=1.2$ nS, $\bar{g}_{\rm syn(inh)\underline{~}max}=0.5$ nS), 2:1 coupling (with inhibition: $\bar{g}_{\rm syn(ext)\underline{~}max}=0.8$ nS, $\bar{g}_{\rm syn(inh)\underline{~}max}=1.5$ nS; no inhibition: $\bar{g}_{\rm syn(ext)\underline{~}max}=0.4$ nS, $\bar{g}_{\rm syn(inh)\underline{~}max}=0.0$ nS), and intermittent ($\bar{g}_{\rm syn(ext)\underline{~}max}=2.0$ nS, $\bar{g}_{\rm syn(inh)\underline{~}max}=3.0$ nS; intermittent regular also corresponds to the same values of synaptic strengths)—observed by changing the strengths of excitatory and inhibitory synaptic conductances. The neurons within each neuronal groups are ranked in ascending order of their bursting activity initiation timings for better visualization. (b) Corresponding averaged population activity of NeuronGroup1 and NeuronGroup2 (top and bottom traces, respectively, in each panel) and (c) Activity of a typical neuron belonging to NeuronGroup1 and NeuronGroup2 (top and bottom traces, respectively, in each panel)

The coupling modes were classified visually based on the averaged population activity. However, since the population activity and the individual neuron activity within the neuronal groups are well correlated, the coupling mode classification could be based on the activity of a pair of neurons from NeuronGroup1 and NeuronGroup2 as well (Fig. 2(c)). Nomenclatures used for classifying the different coupling pattern are as follows:

•  The biphasic mode is characterized by activity pattern of NeuronGroup1 neurons resembling Pre-I neurons—bursting just before inspiratory burst of NeuronGroup2 neurons, inhibition during inspiratory phase, and post-inspiratory second burst (Fig. 2(c), Biphasic).

•  The monophasic coupling mode is characterized by activity of NeuronGroup1 neurons comprising two phases—bursting just before inspiratory burst of NeuronGroup2 neurons and inhibition during inspiratory phase (Fig. 2(c), Monophasic). We termed this coupling mode as ‘monophasic’ because it may be thought of as a reduced form of biphasic coupling where the post-inspiratory second burst of NeuronGroup1 neurons is missing. The bursting frequency of both NeuronGroup1 and NeuronGroup2 neurons are high in monophasic coupling mode; note the magnified time scale used for monophasic coupling in Fig. 2.

•  The synchronous coupling mode is characterized by simultaneous bursting activity of both NeuronGroup1 and NeuronGroup2 neurons (Fig. 2(c), Synchronous).

•  The 2:1 coupling mode is the case where NeuronGroup1 and NeuronGroup2 exhibit bursting activity in 2:1 ratio. 2:1 coupling mode may be further subclassified as (1) 2:1 with inhibition—NeuronGroup1 neuron activity is similar to that of biphasic coupling mode except that its post inspiratory burst is significantly delayed (Fig. 2, 2:1 (with inhibition)), and (2) 2:1 without inhibition—NeuronGroup2 neurons synchronize their bursting activity with every alternate bursting activity of NeuronGroup1 neurons (Fig. 2, 2:1 (no inhibition)). In both of these 2:1 coupling modes, excitatory postsynaptic potentials (EPSPs) are observed on membrane trajectories of NeuronGroup2 neurons coincident with the NeuronGroup1 burst during the expiratory phase.

•  Intermittent coupling mode shows mixed features of synchronous and monophasic coupling mode but with irregular variation in burst duration and inter-burst interval (Fig. 2, Intermittent). Frequently, we also obtained a more regular coupling between NeuronGroup1 and NeuronGroup2 for identical values of excitatory and inhibitory synaptic connections, which we termed as Intermittent (regular)—a pattern in which synchronous and monophasic coupling patterns are alternating (Fig. 2, Intermittent (regular)). Since the occurrence domain of Intermittent (regular) overlaps with that of the conventional Intermittent coupling (see Fig. 3, explanation for which is presented in the following paragraph), we decided to treat these two patterns collectively so as to gain some simplicity in the presentation of results.

Fig. 3

(a) Domain of occurrence of coupling modes as a function of excitatory and inhibitory synaptic strengths of the dual oscillator model. For each combination of $\bar{g}_{\rm syn(ext)\underline{~}max}$ and $\bar{g}_{\rm syn(inh)\underline{~}max}$, seven simulations were performed and the frequency of occurrence of various coupling modes are illustrated using the gray color scale indicated on the right. (b) Representative domain of various coupling modes based on maximum occurrence frequency. For simplicity, the occurrence frequency of the two types of intermittent coupling modes were merged together, since their occurrence domain overlaps

Figure 3(a) shows the domain of various coupling modes observed as a function of $\bar{g}_{\rm syn(ext)\underline{~}max}$ and $\bar{g}_{\rm syn(inh)\underline{~}max}$. Seven simulations were done for each combination of $\bar{g}_{\rm syn(ext)\underline{~}max}$ and $\bar{g}_{\rm syn(inh)\underline{~}max}$ and the occurrence frequencies of each coupling mode are also shown. The boundaries of various coupling mode domains in Fig. 3(a) cannot be precisely determined, because the coupling modes changed slightly every time with the assignment of random initial values for the various synaptic conductance and $\bar{g}_{\rm NaP}$ conductance. Figure 3(b) shows the domain of various coupling modes based on their maximum frequency of occurrence. It is remarked that occasionally the coupling patterns obtained in 2:1 coupling domain could not be meaningfully classified. The values of $\bar{g}_{\rm syn(ext)\underline{~}max}$ and $\bar{g}_{\rm syn(inh)\underline{~}max}$ are low in this region and thus the bursting behavior of NeuronGroup1 and NeuronGroup2 in some simulations were largely uncoupled.

3.2 Simulating the effects of decrease in preBötC neuronal excitability

We simulated the effect of decrease in preBötC neuronal excitability by increasing the leak current conductance $\bar{g}_{\rm L}$ (used for the computation of leakage current I _L (Butera et al. 1999a)) for the neurons constituting NeuronGroup2 and studied its effect on the coupling pattern using the dual oscillator model. We studied the case when $\bar{g}_{\rm syn(ext)\underline{~}max} = 0.8$ nS and $\bar{g}_{\rm syn(inh)\underline{~}max} = 4.5$ nS, which consistently produces a biphasic coupling pattern in our simulation (Figs. 3 and 2(a)). When $\bar{g}_{\rm L}$ is increased, the number of cells with pacemaker property in NeuronGroup2 decreases (Fig. 10). When $\bar{g}_{\rm L}$ is increased beyond 3.4 nS, none of the NeuronGroup2 neurons show pacemaker property. In this situation, NeuronGroup1 is the sole rhythm generator, and NeuronGroup2 is the inspiratory pattern generator triggered by NeuronGroup1 due to the excitatory synaptic connection from NeuronGroup1 to NeuronGroup2 ($\bar{g}_{\rm syn(ext)\underline{~}max}=0.8$ nS).We observed 2:1 coupling, 3:1 coupling and eventually quantal slowing of bursting rhythms between NeuronGroup1 and NeuronGroup2 (Fig. 4), as the value of $\bar{g}_{\rm L}$ was successively increased. In Fig. 4(c), NeuronGroup2 non-deterministically exhibits bursts at the fourth or fifth phasic excitations from the NeuronGroup1. Figure 5 depicts the non-deterministic variability in the interburst interval of NeuronGroup2 for a long duration simulation of quantal phenomenon; the quantal distribution of the interburst interval is readily apparent. However, the quantal slowing phenomenon was not observed in every simulation; the occurrence rate was about 70%; in the rest of the cases we obtained regular 4:1 or 5:1 coupling between NeuronGroup1 and NeuronGroup2, similar to the regular 3:1 coupling depicted in Fig. 4(b). This was because, for every simulation, we assigned different random values to $\bar{g}_{\rm NaP}$, $\bar{g}_{\rm syn(int)}$, $\bar{g}_{\rm syn(ext)}$, and $\bar{g}_{\rm syn(inh)}$; the occurrence of quantal slowing was sensitive to the variability in bursting property of constituent neurons and distributions of the synaptic conductances.

Fig. 4

Effect of increase in $\bar{g}_{\rm L}$ value on the coupling mode between NeuronGroup1 and NeuronGroup2. With the successive increase in $\bar{g}_{\rm L}$, the coupling mode which is originally biphasic (Fig. 2(a), biphasic) changes to (a) 2:1 coupling, (b) 3:1 coupling and eventually (c) quantal slowing—NeuronGroup2 exhibits bursting non-deterministically at fourth or fifth phasic excitations from NeuronGroup1. In raster plots, the neurons within each neuronal groups are ranked in ascending order of their bursting activity initiation timings for better visualization. (d) Averaged population activity of NeuronGroup1 (top trace) and NeuronGroup2 (bottom trace) corresponding to the quantal slowing case depicted in (c)

Fig. 5

(a) Plot depicting the non-deterministic variation of interburst interval of NeuronGroup2 with time. The vertical separation between the horizontal grid lines indicates the interburst interval of NeuronGroup1. NeuronGroup2 exhibits bursts, non-deterministically, at the third, fourth, fifth, sixth, seventh and ninth phasic excitatory drives from NeuronGroup1.(b) Population activities of NeuronGroup1 and NeuronGroup2 for the time interval depicted by the thick gray bar in (a). For the above simulation, the values of $\bar{g}_{\rm syn(ext)\underline{~}max}$, $\bar{g}_{\rm syn(inh)\underline{~}max}$ and $\bar{g}_{\rm L}$ were 0.8 nS, 4.5 nS and 3.89 nS, respectively

Figure 6 depicts the interburst interval of NeuronGroup2 as a function of $\bar{g}_{\rm L}$. To make the quantal slowing effect apparent, the interburst interval of NeuronGroup2 is non-dimensionalized with the interburst interval of NeuronGroup1. There is some variability in the interburst interval of the NeuronGroup1 due to their random initialization. Since, at elevated values of $\bar{g}_{\rm L}$ for neurons in NeuronGroup2, only NeuronGroup1 is rhythm generator and NeuronGroup2 bursts under the influence of excitatory phasic drives from NeuronGroup1, the variability in the interburst interval of NeuronGroup1 is translated to variability in the interburst interval of NeuronGroup2. Non-dimensionlization eliminates this variability dependence and makes the quantal slowing effect pronounced.

Fig. 6

Interburst interval of NeuronGroup2 as a function of $\bar{g}_{\rm L}$. The interburst interval of NeuronGroup2 is non-dimesionlized to remove the variability in interburst interval of NeurnGroup2 caused by variability in the interburst interval of NeuronGroup1. A new simulation is performed for each value of $\bar{g}_{\rm L}$. For $\bar{g}_{\rm L}<$3.84 nS, we observed single bursting frequency of NeuronGroup2. Quantal slowing was observed in the range 3.84 nS $\leq\bar{g}_{\rm L}\leq$ 4.0 nS; multiple data points corresponding to the same value of $\bar{g}_{\rm L}$ in this range depicts the multiple interburst intervals observed. For example, the two solid squares depicts the the two interburst interval of NeuronGroup2 obtained for the case $\bar{g}_{\rm L} = 3.88$ nS shown in Fig. 4(c)

The decrease in preBötC excitability is also simulated by reducing $\bar{g}_{\rm NaP\underline{~}max}$ value of NeuronGroup2. When pursued this course, setting $\bar{g}_{\rm syn(ext)\underline{~}max} = 0.8$ nS and $\bar{g}_{\rm syn(inh)\underline{~}max} = 4.5$ nS in the dual oscillator model, we observed similar changes in coupling mode between NeuronGroup1 and NeuronGroup2 as depicted in Fig. 4 as the $\bar{g}_{\rm NaP\underline{~}max}$ value was lowered; the coupling pattern successively changed from initial biphasic mode to 2:1 coupling, then to 3:1 coupling which eventually lead to quantal slowing. Figure 7 depicts a typical case of quantal slowing obtained by reducing $\bar{g}_{\rm NaP\underline{~}max}$ value.

Fig. 7

Quantal slowing observed when $\bar{g}_{\rm NaP\underline{~}max}$ of NeuronGroup2 was reduced to 2.09 nS ($\bar{g}_{\rm syn(ext)\underline{~}max}$ = 0.8 nS, $\bar{g}_{\rm syn(inh)\underline{~}max}$ = 4.5 nS). (a) Raster plots depicting the bursting activity of neurons within NeuronGroup1 and NeuronGroup2, and (b) Averaged population activity within NeuronGroup1 and NeuronGroup2 (top and bottom trace, respectively). In raster plot, the neurons within each neuronal groups are ranked in ascending order of their bursting activity initiation timings for better visualization

To gain insight into mechanism causing quantal slowing, further investigations were made. The histograms in Fig. 8 show the number of cells in the NeuronGroup2 simultaneously bursting during the quantal slowing case depicted in Fig. 4(c). Almost all cells in NeuronGroup1 showed synchronized activity and hence their activity is not shown. Note that a small fraction of neurons in the NeuronGroup2, variable in number, exhibited bursting by excitatory phasic drives from NeuronGroup1 in the interval between successive synchronized NeuronGroup2 bursting. Next, we removed the interconnections within NeuronGroup2, and investigated the number of neurons that exhibit bursting by phasic drives from NeuronGroup1 under such condition. To do this, we set $\bar{g}_{\rm syn(ext)\underline{~}max} = 0.8$ nS and $\bar{g}_{\rm syn(inh)\underline{~}max} = 0$ nS, and for all NeuronGroup2 cells, we set $\bar{g}_{\rm syn(int)} = 0$.

Fig. 8

Histograms showing the number of cells in NeuronGroup2 exhibiting simultaneous bursting during the quantal slowing case depicted in Fig. 4(c)

Consequently, the dual oscillator model was modified as follows: (1) NeuronGroup1 remains unchanged, (2) NeuronGroup2 becomes essentially a group of unconnected 81 neurons, and (3) the excitatory synaptic connection from NeuronGroup1 to NeuronGroup2 is present ($\bar{g}_{\rm syn(ext)\underline{~}max} = 0.8$ nS) but the inhibitory synaptic connection from NeuronGroup2 to NeuronGroup1 is absent. Subsequently, we successively increased the value of $\bar{g}_{\rm L}$ for neurons in NeuronGroup2 and studied the time history of neurons exhibiting bursting in NeuronGroup2 quantitatively. Figure 9 depicts the details of one representative study when $\bar{g}_{\rm L} = 3.85$ nS for neurons in NeuronGourp2. For the elevated value of $\bar{g}_{\rm L} = 3.85$ nS, none of the neurons in NeuronGroup2 possess pacemaker property (Fig. 10). However, under the influence of excitatory phasic drives from NeuronGroup1, a total of 19 (out of 81) neurons in NeuronGroup2 exhibit bursting at various instances. Moreover, not all 19 neurons in NeuronGroup2 exhibit simultaneous bursting; the histograms in Fig. 9 show the variation in the number of neurons that exhibit simultaneous bursting at each excitatory phasic drive from NeuronGroup1. A maximum of 12 and a minimum of 4 neurons in NeuronGroup2 exhibit simultaneous bursting, indicating that there is considerable variability. The number of neurons in NeuronGroup2 that exhibit bursting with the influence of excitatory phasic drives from NeuronGroup1 for different values of $\bar{g}_{\rm L}$ is shown in Fig. 10.

Fig. 9

Variability in the number of neurons in NeuronGroup2 exhibiting bursting under the influence of excitatory phasic drive from NeuronGroup1. The neurons in NeuronGroup2 are unconnected ($\bar{g}_{\rm syn(int)} = 0$ nS for all neurons in NeuronGroup2) and $\bar{g}_{\rm syn(ext)\underline{~}max}$ = 0.8 nS and $\bar{g}_{\rm syn(inh)\underline{~}max}$ = 0 nS. In raster plot, the neurons within NeuronGroup1 are ranked in ascending order of their bursting activity initiation timings for better visualization

Fig. 10

Typical decrease in the number of neurons with pacemaker property in NeuronGroup2 as $\bar{g}_{\rm L}$ value is increased

4 Discussion

4.1 Coupling modes of pFRG and preBötC rhythm-generating networks

Our model could produce 1:1 coupling of bursting rhythms between pFRG and preBötC with the characteristic biphasic pre-inspiratory and post-inspiratory firing pattern of pFRG neurons. Post-inspiratory burst of NeuronGroup1 neurons during biphasic coupling, characteristic of Pre-I neurons, is the consequence of recovery from I _NaP inactivation during inspiration-related inhibition of the Pre-I firing in our model (post-hyperpolarization rebound bursting (Butera et al. 1999a)).

Our results indicated that the coupling modes depend on the strengths of synaptic connections between the two networks. When we view the coupling modes as a function of the strengths of excitatory input from NeuronGroup1 to NeuronGroup2 ($\bar{g}_{\rm syn(ext)\underline{~}max}$) and inhibitory input from NeuronGroup2 to NeuronGroup1 ($\bar{g}_{\rm syn(inh)\underline{~}max}$), three major domains are recognized: (1) synchronous coupling domain with a low inhibitory strength, (2) biphasic coupling domain with a high inhibitory strength and low/moderate excitatory strength, and (3) monophasic coupling domain with high excitatory and high inhibitory strengths (Fig. 3). The changes in coupling mode by the inhibitory synaptic strength are consistent with experimental results. When a GABA antagonist (bicuculline or picrotoxin) or a glycine antagonist (strychnine) is given to the perfusate of brainstem spinal cord preparations, IIPI becomes absent or negligible, and the activity of Pre-I neurons overlaps with the C4VR inspiratory activity (Onimaru et al. 1990). These findings are in agreement with the model prediction when the strength of inhibitory connection from NeuronGroup2 to NeuronGroup1 is reduced.

Between the synchronous and biphasic major domains, there exists the domain of 2:1 coupling ‘with inhibition’ for moderate values of $\bar{g}_{\rm syn(inh)\underline{~}max}$ (Fig. 3). The 2:1 coupling ‘with inhibition’ may be viewed as a precursor to biphasic coupling as the coupling mode transits from synchronous to biphasic with an increase in inhibitory strength. In the case of 2:1 coupling ‘with inhibition’, the inhibitory strength is not strong enough to cause an immediate post hyperpolarization rebound bursting of NeuronGroup2 as in biphasic coupling; the bursting of NeuronGroup2 is considerably delayed. We observed EPSPs in NeuronGroup2 neuronal membrane trajectories coincident with the rebound bursts of NeuronGroup1. These EPSPs could not evoke NeuronGroup2 burst because of two reasons: (1) low value of $\bar{g}_{\rm syn(ext)\underline{~}max}$, and (2) the timing of excitatory inputs to NeuronGroup2 from its previous burst is too close. With an increase in the value of $\bar{g}_{\rm syn(ext)\underline{~}max}$, excitatory inputs to the neurons of NeuronGroup2 increase, causing it to burst every time the pFRG bursts, thus transforming the 2:1 coupling to the intermittent coupling mode (Fig. 3).

Likewise, the 2:1 coupling ‘without inhibition’, which is primarily found to the left of the synchronous coupling domain and for extremely low values of inhibitory strength, may be viewed as a precursor to synchronous coupling. In 2:1 coupling ‘without inhibition’, the excitatory strength is not strong enough to cause a synchronized bursting of NeuronGroup2 at each phasic excitation from NeuronGroup1; instead it synchronizes at every alternate phasic excitation from NeuronGroup1. It may be noted that the cases of 2:1 coupling ‘with’ and ‘without inhibition’ presented in Fig. 2 lie at the two extremities of the 2:1 coupling domain marked in Fig. 3(b); thus, the difference between them is readily evident. However, the transition from ‘without inhibition’ form of 2:1 coupling to ‘with inhibition’ form is extremely smooth; therefore, we merged the two cases together while determining the coupling domains in Fig. 3.

As 2:1 ‘with inhibition’ coupling mode transits to intermittent coupling mode with increase in the value of $\bar{g}_{\rm syn(ext)\underline{~}max}$, biphasic mode transits to monophasic coupling mode with increase in the value of $\bar{g}_{\rm syn(ext)\underline{~}max}$. For the biphasic coupling mode, the value of $\bar{g}_{\rm syn(ext)\underline{~}max}$ is relatively low. Consequently, the EPSPs provided by NeuronGroup1 at the time of its post inspiratory burst to NeuronGroup2 is not sufficient to make it burst again. However, for higher values of $\bar{g}_{\rm syn(ext)\underline{~}max}$, as in the case of monophasic coupling, the post-inspiratory burst of NeuronGroup1 provides high enough EPSPs to NeuronGroup2, causing it to burst again (Fig. 2, Monophasic). Since the interval between pre-inspiratory and post-inspiratory bursts of NeuronGroup1 is very small, the frequency of bursts in monophasic coupling mode is high. It may be observed that the excitation-inhibition cycle (excitation of NeuronGroup2 by NeuronGroup1 and inhibition of NeuronGroup1 by NeuronGroup2) continues infinitely in monophasic coupling, whereas in case of intermittent coupling, the excitation-inhibition cycle intermittently breaks and restarts again. This suggests that the excitation-inhibition cycle is more stable when both $\bar{g}_{\rm syn(ext)\underline{~}max}$ and $\bar{g}_{\rm syn(inh)\underline{~}max}$ are high.

To illustrate the predictive capability of dual oscillator model, it was remarked earlier in this section that biphasic and synchronous coupling modes are experimentally observed cases. In the same regard, it may also be of significance to note that the 2:1 coupling ‘with inhibition’ as depicted in Fig. 2 and the 2:1 coupling depicted in Fig. 4(a) are also observed in experimental conditions (Okada et al. 2007), which adds to the predictive capability of the model.

4.2 Mechanisms of quantal slowing

During quantal slowing, subthreshold phasic drives to preBötC inspiratory neurons are observed coincident with the timing of skipped inspiratory burst. Based on this observation together with non-deterministic jumps of the bursting period to integer multiples of the control period, Mellen et al. (2003) have suggested that opioid-induced quantal slowing results from transmission failure from unaffected Pre-I neurons to depressed preBötC networks. In addition, they have also mentioned that noisy mutual coupling between the rhythmically active preBötC and the pFRG networks could be an alternative explanation of quantal slowing. Following this, Wittmeier et al. (2008) simulated the quantal slowing phenomenon by incorporating stochastic excitatory synaptic transmission from pFRG to preBötC and reducing its (mean) conductance by 12.5% of control. Thereby, they reinforced the suggestion that quantal slowing results from both transmission failure and noisy mutual coupling between pFRG and preBötC.

In the present study, however, quantal slowing was observed without assuming transmission failure or noisy interactions. Here we use the term ‘transmission failure’ specifically as a suppression of excitatory synaptic transmission from pFRG to preBötC. The results of the present study are both new and simple; ‘simple’ in the sense that the basic dual oscillator model by itself is adequate to simulate quantal slowing without incorporating any additional hypothesis like transmission failure or noisy mutual coupling. Nevertheless, under some experimental conditions, these processes may contribute secondarily to the primary mechanism of quantal slowing as explained in the following paragraph.

The quantal slowing phenomena stems from the fact that the synchronized bursting of NeuronGroup2 is governed by the states of a small fraction of neurons in it that exhibit bursting by phasic excitatory drives from NeuronGroup1 under elevated values of $\bar{g}_{\rm L}$ or decreased values of $\bar{g}_{\rm NaP\underline{~}max}$ for neurons in NeuronGroup2 (Fig. 10). These small fraction of neurons provide excitatory inputs to the other quiescent neurons of NeuronGroup2 and help them to burst. If the total number of neurons in NeuronGroup2 exhibiting bursting exceeds a ‘critical value’, it results in a self-sustained chain of events where progressively larger number of neurons of NeuronGroup2 exhibit bursting and eventually all the neurons in NeuronGroup2 exhibit synchronized bursting. On the other hand if the number of bursting neurons in NeuronGroup2 falls short of the critical value, then the neurons in NeuronGroup2 do not exhibit synchronized bursting. Figure 9 depicted variability in the number of neurons in NeuronGroup2 exhibiting bursting by excitatory phasic drives from NeuronGroup1. Though, the case presented in Fig. 9 corresponds to unconnected neurons in NeuronGroup2, it is likely that similar variability, albeit to a different degree, exists even when the neurons in NeuronGroup2 are interconnected through chemical synapses. This is evident if one compares the similarity in the bursting activity of a small fraction of NeuronGroup2 neurons in the interval between synchronized bursting in Fig. 8 (and Fig. 7) with that of the uncoupled neurons in NeuronGroup2 depicted in Fig. 9. Consequently, the number of bursting neurons in NeuronGroup2 exceeds the critical number intermittently and synchronized NeuronGroup2 bursting occurs. This results in quantal slowing of synchronized bursting of NeuronGroup2. From our simulation study, we observed that a maximum of about 15 neurons in NeuronGroup2 exhibit bursting when NeuronGroup2 fails to exhibit synchronized bursting. Thus, the ‘critical number’ of bursting neurons that necessarily be surpassed for synchronized NeuronGroup2 bursting may be greater than 15 (out of 81) neurons. The critical number is expected to be fuzzy rather than precisely defined. The fuzziness stems from the fact that the distribution of $\bar{g}_{\rm NaP}$ and $\bar{g}_{\rm syn(int)}$ among NeuronGroup2 neurons is randomly assigned. Moreover, the critical number is a function of the network topology as well (see Section 4.3). For a neuronal group where the constituent neurons are interconnected in a different way than what we have considered here, the critical number is different. However, the qualitative aspects of the results, such as intermittent failure of NeuronGroup2 to exhibit synchronized bursting, are generally expected and could be an underlying mechanism of quantal slowing.

Quantal slowing was not observed in 30% of the simulations. In these cases, the distributions of random values of $\bar{g}_{\rm NaP}$ and $\bar{g}_{\rm syn(int)}$ within NeuronGroup2 are such that its constituent neurons are relatively better synchronized. This was evident from visual observation of the raster plots obtained for these cases, which resembled the regular coupling pattern in Fig. 4(b). Consequently, there is no intermittent failure of NeuronGroup2 to exhibit synchronized bursting and hence no quantal slowing is observed. In such cases, the activity of the neuronal groups, both NeuronGroup1 and NeuronGroup2, may be essentially captured by replacing each of them with a single pacemaker neuron. This probably explains why Wittmeier et al. (2008) could not simulate quantal slowing phenomenon with their simplified model unless they incorporated stochastic synaptic transmission. The simplified model lacks the complexity necessary to capture the essential source of (apparent) non-determinism resulting the quantal slowing phenomenon (as discussed in the previous paragraph), which is intrinsic in the neuronal population model that we have considered. It is important to realize that the non-determinism associated with the occurrence of quantal slowing simulation is an ‘apparent’ one—the dual oscillator model used to simulate quantal slowing is deterministic; it is the nonlinearity present in the bursting behavior of individual neurons within the neuronal groups that result in non-deterministic-like intermittent synchronized bursts of NeuronGroup2. Thus, our results also hint toward the fact that quantal slowing may not necessarily be derived from any stochastic processes but may actually be a deterministic phenomena.

Figures 5 and 6 suggests that the interburst durations of NeuronGroup2 while exhibiting quantal slowing are not integer multiples of the control period but are instead fractional multiples. Interestingly, this was elucidated by Wittmeier et al. (2008) as well, though they achieved it by incorporating a stochastic synaptic transmission. Consequently, it may be inferred that the fractional nature of quantal slowing is primarily a characteristic of dual oscillator models and is independent of the mechanism inducing the non-determinism.

4.3 Relation of preBötC neuronal network topology to quantal slowing

We have assumed that each neuron of NeuronGroup2 is synaptically connected to every other neuron within NeuronGroup2 (“all-to-all” connectivity). Consequently, the intra-neuronal group synaptic current I _syn(int) for individual neurons of NeuronGroup2 (the first term on the right hand side of Eq. (2)) is:

$$ I_{\rm syn(int)}=\left(\sum\limits_{l\in L}g^{l}_{\rm syn(int)}\cdot s^{l}_{\rm int}\right)\left(V-E_{\rm syn(e)}\right) $$

(8)

where L is NeuronGroup2. However, we subsequently reduced the above equation to

$$ I_{\rm syn(int)}=\bar{g}_{\rm syn(int)}\cdot\bar{s}_{\rm int}\left(V-E_{\rm syn(e)}\right) $$

(9)

(refer Section 2.2, in particular the reduction of Eq. (2) to Eq. (3)) and then randomly assigned values to $\bar{g}_{\rm syn(int)}$ from a uniformly distributed probability density function ranging between 0 and 10 nS. Since, $\bar{s}_{\rm int}$ is the same for all the neurons, it is primarily the variability in the values of $\bar{g}_{\rm syn(int)}$ across the neurons in NeuronGroup2, which contributes to the variability in the intra-neuronal synaptic currents for individual neurons of NeuronGroup2. From synchronization perspective, a neuron with low value of $\bar{g}_{\rm syn(int)}$ may be considered as weakly connected to rest of the neurons within NeuronGroup2 and vice versa.

It is remarked that the reduction of Eq. (8) to Eq. (9) modifies the intrinsic homogeneity present in “all-to-all” network connectivity. If we randomly assign values to $g^{l}_{\rm syn(int)}$ from a uniformly distributed probability density function of any appropriate range, and perform the summation shown in Eq. (8), the variability in the individual values of $g^{l}_{\rm syn(int)}$ will be equalized resulting in almost homogeneous distribution of I _syn(int) across the neurons of NeuronGroup2. Homogenous distribution of I _syn(int) within NeuronGroup2 would imply that each neuron is connected equally strongly to rest of the neurons within NeuronGroup2. Since this is unlikely to be the case in the real preBötC, we presume that our way of introducing variability in the distribution of I _syn(int) across the neurons of NeuronGroup2 is qualitatively more realistic. The significance of the usage of uniformly distributed probability density functions for randomly assigning the values to $\bar{g}_{\rm syn(ext)}$ and $\bar{g}_{\rm syn(inh)}$ in the model may be interpreted similarly.

Now, since in our simulations, quantal slowing is observed primarily due to intermittent failure of NeuronGroup2 to exhibit synchronized burst, it is remarked that the ability of NeuronGroup2 to exhibit synchronized bursts is intrinsically dependent on its network topology, in particular, the distribution of I _syn(int) across the neurons of NeuronGroup2. In the quantal slowing simulation, Fig. 4(c) for example, we observe that NeuronGroup2 mostly exhibits synchronized bursts at every fourth phasic excitation from NeuronGroup1 but frequently it skips the fourth phasic excitation from NeuronGroup1 and bursts at the fifth phasic excitation. Thus, once NeuronGroup2 has exhibited a synchronized burst, there is an element of (apparent) non-determinism whether it will exhibit its next synchronized bust at the fourth phasic excitation from Neurongroup1 or not. However, in experimentally observed quantal slowing, preBötC elucidates non-determinism in its busting at each phasic excitation of pFRG (Mellen et al. 2003). We reason that this discrepancy is probably due to the difference in the neuronal network topology of real preBötC as compared to the one we have incorporated in NeuronGroup2.

Thus, it is remarked that the highlight of the present simulation study of quantal slowing is only its essential qualitative feature—the (apparent) non-determinism associated with bursting of preBötC—and the potential mechanism causing it. The quantitative features of quantal slowing may only be simulated by incorporating a more realistic network topology for NeuronGroup2.

4.4 Some additional simulation studies

In the model presented above, we introduced variability in bursting frequencies of the constituent neurons of NeuronGroup1 and NeuronGroup2 by providing variability in the allocated value of persistent Na⁺ current conductance $\bar{g}_{\rm NaP}$ of neurons. The variability in bursting frequencies may also be introduced by providing variability in the leakage current conductance $\bar{g}_{\rm L}$ of neurons keeping the $\bar{g}_{\rm NaP}$ value fixed. We pursued this course too and verified that all the qualitative features presented above remains unchanged. For this study, we randomly assigned $\bar{g}_{\rm L}$ to neurons of each neuronal group from a uniform probability distribution function ranging from 2.3 to 2.8 nS; we set $\bar{g}_{\rm NaP} = 3.5$ nS for all neurons of NeuronGroup1 and $\bar{g}_{\rm NaP} = 2.5$ nS for that of NeuronGroup2, so that the bursting frequency of NeuronGroup1 is higher than that of NeuronGroup2.

Further, we used a Gaussian distribution (with a standard deviation of 0.3 nS about the mean values of $\bar{g}_{\rm NaP}=3.5$ nS for NeuronGroup1 and $\bar{g}_{\rm NaP}=2.45$ nS for NeuronGroup2) instead of the uniform distribution, for allocating the values of $\bar{g}_{\rm NaP}$ to individual neurons. Again the qualitative features of results obtained were essentially the same as presented above. Thus, it is suggested that the results presented in the manuscript are independent of the mode of introducing the variability in bursting frequencies of the constituent neurons of the neuronal groups.

5 Model justification

We assumed that the intrinsic bursting frequency of pFRG is faster than that of preBötC. This is a prerequisite for the consistent firing of pFRG neurons before the burst of preBötC neurons. To our knowledge, this assumption has never been tested experimentally. The intrinsic bursting rhythm of pFRG neurons, separated from preBötC can be monitored from the facial nerve rootlet, however, the inherent bursting frequency of the pFRG is difficult to determine because the rhythm is suppressed by pons or other brain regions with the facial nerve attached (Onimaru et al. 2006). In the dual oscillator model, we simulated by reversing the bursting frequencies of NeuronGroup1 and NeuronGroup2, that is by setting $\bar{g}_{\rm NaP\underline{~}max} = 3.0$ nS for NeuronGroup1 and $\bar{g}_{\rm NaP\underline{~}max} = 4.0$ nS for NeuronGroup2. We did not obtain the characteristic biphasic coupling for any combination of excitatory and inhibitory synaptic strengths (the range of synaptic strengths considered for this study was the same as that in Fig. 3).

We assumed that the persistent Na⁺ current is responsible for rhythmic bursts of pFRG and preBötC neurons. This assumption has been recently contradicted. Pace et al. (2007a) have demonstrated that I _NaP does not contribute to inspiratory drive potential generation in the vast majority of preBötC neurons. Another intrinsically bursting neuronal group has been identified in the preBötC, which is dependent on a calcium-activated nonspecific cationic current (I _CAN) (Peña et al. 2004). However, the relative contribution of pacemaker properties and synaptic inputs to rhythm generation is likely to be dynamic (Johnson 2007). For example, during hypoxia, the respiratory pattern generator is more dependent on pacemaker properties and produces gasping behavior (Paton et al. 2006). Further recent data of Koizumi et al. (2008) from neonatal rats contradict the data of Pace et al. (2007b) from mice, and show that inspiratory rhythm generation can be I _NaP-dependent in the isolated preBötC network in vitro. Therefore, although we did not consider I _CAN, the breakdown of synchronized burst could still be the underlying mechanism of quantal slowing in certain conditions.

It has been shown that glutamatergic excitatory synaptic inputs are required to evoke the I _CAN-dependent inspiratory drive potentials (Pace et al. 2007b), suggesting that the rhythmic burst is generated by a “group pacemaker” mechanism (Feldman and Del Negro 2006) in which rhythm generation is an emergent property of the network. Recently, Rubin et al. (2009) proposed a mathematical model of the neurons exhibiting “group pacemaker” mechanism. However, this model lacks the post-hyperpolarization rebound bursting dynamics which is essential for producing the post-inspiratory burst of pFRG, and thus we cannot use this neuronal model to simulate and study quantal slowing.

We simulated the opioid-induced reduction of excitability in preBötC bursting neurons both by increasing $\bar{g}_{\rm L}$ and by reducing $\bar{g}_{\rm NaP\underline{~}max}$. However, we do not know how opioids depress preBötC neurons. Opioids activate a G-protein-coupled inwardly rectifying potassium conductance known as G _irk , resulting in the hyperpolarization of neurons throughout the central nervous system (Williams et al. 1982; Wimpey and Chavkin 1991). Therefore, the depressant effect of opioids on preBötC may be reasonably simulated by an increase in $\bar{g}_{\rm L}$. On the other hand, it is unlikely that opioids directly affect $\bar{g}_{\rm NaP\underline{~}max}$ because pFRG bursting neurons, which are I _NaP-dependent (Onimaru et al. 1997), are not affected by opioids (Takeda et al. 2001; Janczewski et al. 2002).

Since the two putative rhythm-generating networks are embedded hierarchically in the central respiratory pattern-generating network in more intact animals (Smith et al. 2007), dynamic interactions among neuronal groups must not be so simple. Opioid-induced quantal slowing has been observed in juvenile (Mellen et al. 2003) and adult (Vasilakos et al. 2005) rats in vivo. These observations can be interpreted as the manifestation of dynamic interactions between fundamental-level networks. However, it has been reported that eupnea of in situ intra-arterially perfused rats persists following the blockade of the two burst-generating currents, I _NaP and I _CAN (St.-John 2008). Therefore, whether the dual oscillator configuration is the fundamental-level network component of eupnea or that of gasping remains to be clarified.

In summary, we developed a dual oscillator model to understand the dynamic interactions between pFRG and preBötC neurons. Our model essentially assumes (1) both pFRG and preBötC networks are rhythm generators, (2) preBötC receives excitatory inputs from pFRG, and pFRG receives inhibitory inputs from preBötC, and (3) persistent Na⁺ current conductance and synaptic current conductances are randomly distributed. Our model could produce the characteristic behaviors observed experimentally in neonatal brainstem spinal cord preparations. The coupling mode depended on the strengths of excitatory and inhibitory connections of the oscillator. In contrast to the earlier suggestions, quantal slowing was observed without transmission failure (suppressed excitatory synaptic conductance) or noisy mutual interactions between the neuronal networks of the oscillator. Our study suggests that quantal slowing may actually be a deterministic phenomena and the non-determinism associated with it may only be an ‘apparent’ one. We suggest that quantal slowing results from inhomogeneous properties of individual cells within the oscillator and subsequent breakdown of synchronized bursting within the preBötC oscillator.