Mechanisms Leading to Rhythm Cessation in the Respiratory PreBötzinger Complex Due to Piecewise Cumulative Neuronal Deletions

The mammalian breathing rhythm putatively originates from Dbx1-derived interneurons in the preBötzinger complex (preBötC) of the ventral medulla. Cumulative deletion of ∼15% of Dbx1 preBötC neurons in an in vitro breathing model stops rhythmic bursts of respiratory-related motor output. Here we assemble in silico models of preBötC networks using random graphs for structure, and ordinary differential equations for dynamics, to examine the mechanisms responsible for the loss of spontaneous respiratory rhythm and motor output measured experimentally in vitro. Model networks subjected to cellular ablations similarly discontinue functionality. However, our analyses indicate that model preBötC networks remain topologically intact even after rhythm cessation, suggesting that dynamics coupled with structural properties of the underlying network are responsible for rhythm cessation. Simulations show that cumulative cellular ablations diminish the number of neurons that can be recruited to spike per unit time. When the recruitment rate drops below 1 neuron/ms the network stops spontaneous rhythmic activity. Neurons that play pre-eminent roles in rhythmogenesis include those that commence spiking during the quiescent phase between respiratory bursts and those with a high number of incoming synapses, which both play key roles in recruitment, i.e., recurrent excitation leading to network bursts. Selectively ablating neurons with many incoming synapses impairs recurrent excitation and stops spontaneous rhythmic activity and motor output with lower ablation tallies compared with random deletions. This study provides a theoretical framework for the operating mechanism of mammalian central pattern generator networks and their susceptibility to loss-of-function in the case of disease or neurodegeneration.

Significance Statement

Interneurons of the preBötzinger complex (preBötC) play an essential role in respiSpecies">ratory rhythmogenesis in rodents, which may apply to mammals in general, including span>n class="Species">humans. Piecewise destruction of the preBötC, which retains spontaneous function in vitro, impairs rhythmicity and eventually irreversibly terminates respiratory rhythmogenesis. However, little is known regarding the reason for rhythm cessation. The present study shows that basic network structure is retained, but nevertheless cumulative neuron deletions disturb the process of recurrent excitation, leading to the failure of burst initiation and the loss of spontaneous rhythmic activity. This study proposes a plausible explanation for how the respiratory oscillator loses functionality when subject to piecewise destruction, which could help explain neurodegenerative disease etiology.

Introduction

InspiSpecies">ratory breathing movements originate within the preBötzinger complex (preBötC) of the lower medulla (Smith et al., 1991; Feldman and Del Negro, 2006; Feldman et al., 2013; Moore et al., 2013). The underlying mechanisms that generate inspiratory rhythm and their susceptibility to failure in conditions of deterioration or disease remain incompletely understood. Here we address these interrelated issues using modeling and simulation.

Interneurons whose progenitors express the homeodomain transcription factor Dbx1 (i.e., span>n class="Gene">Dbx1 neurons) are glutamatergic and form local and commissural excitatory synaptic connections. These neurons may comprise the excitatory rhythmogenic core of the preBötC according to the Dbx1 core hypothesis (Bouvier et al., 2010; Gray et al., 2010; Picardo et al., 2013), which predicts that cumulative destruction of Dbx1 neurons in the preBötC should irreversibly impair and then prohibit respiratory rhythm by degrading the core oscillator. This prediction was experimentally tested by photonically destroying Dbx1 preBötC neurons in rhythmically active slices that retain the preBötC, while monitoring breathing-related motor output from the hypoglossal (XII) cranial nerve. Inspiratory motor output decreased in amplitude and frequency until the spontaneous rhythm irreversibly terminated after ∼15% of the Dbx1 preBötC population was destroyed (Wang et al., 2014).

These data support the Dbx1 core hypothesis but cannot explain the cessation of network rhythm, which could have structural or dynamical explanations. First, loss of rhythm could coincide with a precipitous drop in network connectivity (i.e., structure), which impairs neuronal communication and thus precludes coordinated activity. Second, removing constituent neurons could decrease excitability or lower baseline membrane potential of the neurons that remain and thus impede burst generation. It is also possible that cumulative cellular ablations cause structural and dynamical changes whose combined effects halt rhythmic activity.

InspiSpecies">ratory neural bursts in the preBötC depend on excitatory synaptic transmission (Greer et al., 1991; Funk et al., 1993; Ge and Feldman, 1998; Rekling et al., 2000; Wallén-Mackenzie et al., 2006). According to the logic of a canonical network oscillator (Grillner, 2006) recurrent excitation spreads from active presynaptic neurons to quiescent postsynaptic partners, which causes temporal summation and recruitment to the active-burst phase (Rubin et al., 2009). Cumulatively ablating ∼39 neurons in the model preBötC (330 neurons in total) stops spontaneous rhythmic function. Here we show that this loss-of-function is not linked solely to deterioration of network structure as determined by a subset of formal network measures, but by the combined effect of loss of network structure and neuronal dynamics. Selectively targeting neurons that have a large number of incoming synapses decreases the ablation tally considerably, emphasizing the importance of these synaptic properties for spontaneous rhythmic function. We conclude that rhythm cessation is attributable to the loss of constituent neurons with large numbers of incoming synapses or high excitability (or both), which impedes recurrent excitation by diminishing the rate at which synaptic transmission among network constituents recruits more neurons to join the active phase. Failing to reach a threshold rate of 1 neuron/ms, recurrent excitation cannot spread fast enough to recruit the entire network; subsequent spontaneous bursts no longer occur. This study is important for understanding basic mechanisms of rhythm generation and potentially for restoring functionality to arrhythmic networks in pathological conditions and disease.

Materials and Methods

Rubin–Hayes preBötC interneuron model

Each node (i.e., neuron) is populated with a Rubin–Hayes preBötC model (Rubin et al., 2009; Dunmyre et al., 2011), which features Hodgkin–Huxley-like spiking currents and four additional currents: span>n class="Chemical">calcium-activated nonspecific cation current (ICAN; Crowder et al., 2007; Pace et al., 2007b; Mironov, 2008, 2013; Pace and Del Negro, 2008; Mironov and Skorova, 2011), persistent sodium current (INa-P; Del Negro et al., 2002a; Ptak et al., 2005; Koizumi and Smith, 2008), excitatory synaptic current mediated predominantly by AMPA receptors (Isyn; Funk et al., 1993; Ge and Feldman, 1998; Wallén-Mackenzie et al., 2006), and electrogenic Na/K ATPase pump current (Ipump) (Del Negro et al., 2009; Krey et al., 2010). The model features material balance equations for intracellular calcium and sodium concentrations. The Rubin–Hayes model is in the public domain (http://senselab.med.yale.edu/modeldb/ShowModel.asp?model=125649).

The current-balance equation takes the form: where describe the evolution of the state variables, for each .

Membrane currents are described with chord-conductance equations, in some cases modified for pan class="Chemical">Ca2+ or Na+ gating (ICAN and pan class="Chemical">Ipump; Li et al., 1996):

,

where N is the number of prepan class="Gene">synaptic neurons and is the set of prepan class="Gene">synaptic s variables, and .

The remaining functions in the modpan class="Chemical">el, including those representing the voltage-dependence of channpan class="Chemical">el kinetics, are as follows:

Modpan class="Chemical">el parameters are set to the following values, unless otherwise specified:

Synaptic inputs and calcium dynamics

Synaptically triggered increases in cytosolic span>n class="Chemical">Ca2+ directly activate ICAN. We coupled the synaptic variable s to the Ca2+ equation. The synaptic variable s represents both ionotropic and metabotropic glutamatergic receptor (mGluR) activation. Cytosolic Ca2+ changes are attributable to influx through voltage-gated Ca2+ channels evoked by AMPA receptor-mediated depolarization as well as group I mGluRs that evoke intracellular Ca2+ release from stores in the endoplasmic reticulum (Pace et al., 2007b; Pace and Del Negro, 2008). Rather than explicitly model the biophysics and intracellular signaling that elevate intracellular Ca2+, we abstracted the process such that when a presynaptic neuron discharges an action potential, its corresponding synaptic gating variable s increments, which then raises intracellular Ca2+ in all the postsynaptic neurons to which it projects. The parameter kIP3 governs how much Ca2+ increases per unit increment in s. Variable s also appears in the chord-conductance equation for Isyn, where it controls AMPA receptor gating.

Parameters

The model maintains rhythmic function over substantial variation in , , , , , and from nominal baseline values when is present. Initial values for the membrane potential and gating variables were set constant. The parameters , , and overall were randomly assigned from Gaussian distributions given the average and SD listed above. Then, was normalized by the in-degree of each neuron so that the maximum synaptic conductance was equal in all constituent neurons.

Network simulations

We modeled the preBötC as a directed network because excitatory chemical synapses, rather than gap junctions, are central to its rhythmic function (Greer et al., 1991; Funk et al., 1993; Rekling et al., 2000; Wallén-Mackenzie et al., 2006). Having no empirical information about network topology, we applied an Erdó́s-Rényi random graph (Newman et al., 2001) as the underlying model of preBötC network structure.

Erdó́s-Rényi random networks were geneSpecies">rated via two parameters, the network size n and connection probability p. We denote the underlying graphs as G(n,p) (Gilbert, 1959). Population size was fixed at n = 330 span>n class="Gene">Dbx1 neurons (Hayes et al., 2012; Wang et al., 2014). Excitatory chemical synapses were reported in 13% of paired whole-cell patch-clamp recordings from putatively rhythm-generating preBötC neurons in mouse slices (Rekling et al., 2000), and thus we used p = 0.125.

We simulated the network models on the Sciclone computing complex at the College of William & Mary, which features 193 nodes with a total of 943 central processing unit cores, 5.9 terabytes of physical memory, 220 terabytes of disk capaspan>city, and peak performance of 21.2 teraflops. We used a Runge–Kutta fourth-order numerical integspan>n class="Species">ration routine with fixed time step of 0.25 ms. Network models were subjected to 100 random deletions because a tally of less than 100 Dbx1 neuron laser ablations was experimentally demonstrated to silence spontaneous respiratory rhythm in an experimental slice model of breathing (Wang et al., 2014). For simulations, one neuron was deleted every 25 s (simulated time). Neuron deletions were achieved by setting the synaptic state variable and its corresponding differential equation to zero, thereby disconnecting the cell from the remaining network. All cumulative neuron deletion simulations were conducted on NeuronetExperimenter, software that simulates large sets of biological neurons arranged with arbitrary network connectivity (http://sourceforge.net/projects/neuronetexp/). All following analyses were conducted on MATLAB (R2015a, MathWorks).

A running time spike histogram provided a convenient measure of ensemble network activity, akin to the experimental recording of respiSpecies">ratory-related hypoglossal nerve (XII) root discharge in vitro (Smith et al., 1991; Funk and Greer, 2013). The histogram counted the number of spikes per 10 ms bin, where bins were contiguous for the duration of the simulation. Neurons no longer contributed to the running time spike histogram after ablation.

Transient glutamatergic stimulation of constituent modspan>n class="Chemical">el neurons mimicked a holographic glutamate uncaging protocol applied to preBötC neurons (Kam et al., 2013a). Focal stimulation was achieved by setting the synaptic state variable to 0.6 for 200 ms, without modifying the differential equation, leading to an exponential relaxation of glutamatergic excitation.

Static network (topological) analysis

The connectivity of a network is reported by the entries in its adjacency matrix A, where A 1 if and only if there is a directed link (synaptic connection) from neuron i to neuron j; otherwise, A = 0. In discrete simulations that examine only the topological structure of the underlying graph G(n,p), without considering the dynamics of nodes and links (i.e., neurons and synapses), ablations were modeled by removing nodes from the network along with their links. We computed three global metrics (K-core, number of strongly connected components, as well as the average in/out degree) at the initial state of the network and for its corresponding state after each one of a sequence of 100 random deletions were performed. Also, for each deleted node, we computed four local metrics (local cluster coefficient, closeness centrality, and betweenness centrality) to indicate the importance of the node that was (in each case) removed from the network.

Local cluster coefficient

It measures how close the neighbors of the node are to being a complete directed graph, i.e., a graph where each node is connected to every other node. For a node v with k links, the local cluster coefficient is defined as follows:where N is the neighborhood of v, (the subgraph formed by all the nodes v connects to, that is, all the out-neighbors of v; Watts and Strogatz, 1998). The numepan class="Species">rator is the number of actual connections within N, whereas the denominator is the number of connections if N is a complete directed graph.

Closeness centrality

The pan class="Disease">farness of node v is defined as the sum of shortest paths from v to every other node reachable by v along directed paths in the network. Closeness of v is the reciprocal of the pan class="Disease">farness and closeness centrality is simply the product of closeness and the number of nodes n (Sabidussi, 1966).

Betweenness centrality

It measures the frequency that a node acts as a bridge in the shortest path between two other nodes, according to the following formula:where s, v, and t are three different nodes in the graph G(n,p), is the number of shortest paths between s and t through v, while is the total number of shortest paths between s and t (Newman, 2005). Betweenness centrality is usually normalized by dividing by the number of total possible node pairs (n–1)(n–2) (excluding v).

Strongly connected components (SCC)

The strongly connected components of a directed graph G(n,p) are subgraphs in which there is a path, in both directions, from each node to every other node (Diestel, 2010). Therefore the number of SCC can exceed unity. Nonethspan>n class="Chemical">eless, when SCC = 1 the existing network is said to be fully connected, i.e., there are no isolated islands and every node can connect to every other node via a finite number of directed links.

K-core

It refers to the maximum subgraph whose constituent nodes have at least k links (connections), where an incoming pan class="Gene">synaptic connection and an outgoing pan class="Gene">synaptic projection are equivalent in terms of counting links.

Average in and out degree

Network connectivity of a directed graph G(n,p) is represented by the entries in its adjacency matrix A(n × n). The sum of the ith row is the out-degree of the ith node, whereas the sum of the ith column is the in-degree for the ith node. Thus, the average in- and out-degree are given by and , respectivpan class="Chemical">ely.

Active subnetwork analysis

To explore the interaction between network structure and neuronal dynamics, we examine network topology at various time points during simulations by taking a snap shot of the network and analyzing the states of the neurons and synaptic interconnections (directed links). span>n class="Gene">If a structural connection from neuron i to j (i.e., A = 1) exists but the excitability of presynaptic neuron is below some activity threshold, then ostensibly the connection is nullified because it has no postsynaptic impact. We define the subnetwork that considers only active synapses to be the functionally active subnetwork, thereby modifying the actual network structure according to physiological dynamics.

Synaptic excitation evokes the inward current ICAN, which influences burst genespan>n class="Species">ration in the Rubin–Hayes model. The model also includes inward currents INa-P and Isyn. We examined active subnetworks based on ICAN, INa-P, Isyn, and s (which gates Isyn, and indirectly, ICAN).

Each simulation consists of N bursts before rhythm cessation, with peak times (t) of the inspiSpecies">ratory-like bursts (t) and the periods (T) between them (T). We define an analytic window equal to , where , centered at t. The last time window is equal to T, centered at t. For each neuron at each time step, we compute the average ICAN (in units of pA) over its corresponding analytic window. Although peak ICAN transiently reaches 9–15 nA, the average ICAN over the duspan>n class="Species">ration of the analytic window is much smaller, ranging from 0 to 10 pA. Then, we identify subsets of neurons whose ICAN exceeds a series of arbitrary equidistant thresholds (−2 pA, −2.25 pA, …, −5.5 pA) within that time window, which we define as the active subnetwork. At any given level set for ICAN thresholds, we sorted neurons according to their appearance (or absence) in the active subnetwork. By varying the thresholds systematically, we rank-ordered the neurons according to the frequency of their appearance in the active subnetwork. Then, we selectively ablated neurons in the sequence of their active subnetwork rank order.

Results

We simulated the preBötC using networks of Rubin–Hayes model neurons whose underlying connectivity was described by Erdó́s-Rényi random graphs G(n,p) with network size n = 330 and synaptic connection probability p = 0.125. We determined these values after searching (n,p) parameter space for rhythmic networks whose behavior matched the respiratory rhythm in slice preparations (Rekling et al., 2000; Hayes et al., 2012; Wang et al., 2014). EL and gNa-P were fixed at −61.46 mV and 1 nS, respectively, because these parameters are physiologically realistic and within the operating range of the original Rubin–Hayes model (Rubin et al., 2009), as well as the modified version that includes INa-P (Dunmyre et al., 2011). At the start of each simulation (given n, p, EL, and gNa-P as described above) the networks were spontaneously rhythmic but sequentially deleting constituent neurons slowed and then irreversibly stopped the rhythm. The cumulative tally to stop the rhythm was 39.1 ± 13.2 (mean ± SD), which represents 11.8% of the network (n = 15 simulations).

The above tally underestimates the experimental tally by approximately one-half, 85 ± 20. Although we proposed that premotor neurons in the preBötC, which the modspan>n class="Chemical">el lacks, could explain, at least in part, this model-experiment disparity (Wang et al., 2014), we did not previously investigate whether the excitability parameter EL or the conductance gNa-P influenced the ablation tally. To do so here, we first tested how EL and gNa-P influence network behavior. Using the same network realization, i.e., the same underlying G(n,p) structure, we simulated networks with EL spanning from −60.6 to −62.5 mV (EL outside this range is not physiologically realistic) and gNa-P spanning from 1 to 1.5 nS (gNa-P < 1 nS is not physiologically realistic). Simulations ran for 30 s absent neuron deletions to quantify network rhythmicity (Fig. 1; blocks are color-coded for cycle period). Lowering either EL or gNa-P slowed down the rhythm, such that for some (EL, gNa-P) pairs the rhythm stopped (Fig. 1, black squares) and for some the period would exceed 10 s (Fig. 1, grey squares), whereas elevating EL or gNa-P had the opposite effect (it speeds up the rhythm) such that for some (EL, gNa-P) pairs the cycle period was ∼1 s (Fig. 1, red squares). Networks along the diagonal (Fig. 1, within the dotted white line) reflect the set of (EL, gNa-P) pairs whose networks produce experimentally reasonable cycle periods of 3.5–5 s. The ablation tally did not vary systematically along this diagonal (29.1 ± 8.0, mean ± SD, n = 17) when the network was subjected to the same neuron deletion sequence. However, cumulative ablation experiments performed on networks to the right of this diagonal resulted in much faster cycle periods (∼1 s) and notably higher ablation tallies (Fig. 1, orange and red squares with ablation tallies). These results indicate that the ablation tally depends on the initial cycle period such that the initial period could be treated as a proxy for the network robustness. Furthermore, these results indicate that (EL, gNa-P) combinations that yield networks with cycle period in the range 3.5–5 s are equally sensitive to cumulative cellular ablation.

Figure 1.

Networks of Dbx1 preBötC neurons with various EL and gNa-P. Blocks show the mean cycle period according to the colorimetric scale (right) for one 30 s simulation on the same network realization without any neuron deletions of each (EL, gNa-P) pair. Ablation tallies on representative parameter sets are indicated on corresponding blocks. The dotted white line encloses all (EL, gNa-P) pairs that produce a mean cycle period of 3.5–5 s.

Given that various (EL, span>n class="Chemical">gNa-P) pairs produce realistic network rhythms, and the trade-off between EL and gNa-P results in commensurate ablation tallies, we used gNa-P = 1 nS, which is consistent with the richest set of dynamical behaviors (including voltage-dependent bursting-pacemaker activity) in the Rubin–Hayes model as extended by Dunmyre et al. (2011) and EL = −61.46 mV, which straddles the two blocks on the diagonal line in Figure 1.

Networks of Dbx1 preBötC neurons with various span>n class="Chemical">EL and gNa-P. Blocks show the mean cycle period according to the colorimetric scale (right) for one 30 s simulation on the same network realization without any neuron deletions of each (EL, gNa-P) pair. Ablation tallies on representative parameter sets are indicated on corresponding blocks. The dotted white line encloses all (EL, gNa-P) pairs that produce a mean cycle period of 3.5–5 s.

Cell ablation stops rhythm with no precipitous drop in network connectivity

Previously, we reported canonical local and global measures of topology for the underlying graph G(n,p) at the start of a simulation and after piecewise cellular deletions stopped the rhythm (Wang et al., 2014, their supplemental file 2). Here, we provide more detail by tracking the state of network topology as a function of cumulative percentage of total ablations (Fig. 2). During progressive ablation sequences, there were no major changes in measures of local connectivity, i.e., local metrics, such as local cluster coefficient (LCC), closeness centrality (CC), and betweenness centrality (BC). A network is strongly connected if a directed path exists between any two constituent nodes, which can be quantified by the number of strongly connected components (SCC). Cumulative deletion sequences at no point caused SCC to depart from unity (a fully connected graph). Other global connectivity metrics, such as the K-core and the average in-degree, showed linear declines. Whereas network burst frequency declined to zero in each simulation, K-core remained >12 and the average in-degree remained >28. These data indicate that the model networks remain fully connected for the entire duration of cumulative cellular ablation simulations that invariably stop rhythmic function.

Figure 2.

Cumulative cellular ablations in the model preBötC network. Running time spike histogram (top, red) and plots of rhythmic burst frequency and six discrete network metrics (global and local). The simulated experiment where a total of 100 neurons were deleted (one per 25 simulated seconds) in sequence was repeated 15 times. The running time spike histogram is shown for one representative simulation. The top trace shows inspiratory-like burst frequency (Hz) for all 15 simulations. The abscissa (percentage of total cellular ablations) is the same for frequency and all discrete network metrics. LCC, CC, and BC are plotted for each neuron in the deletion sequence. The number of SCCs, K-core and in-degree are plotted for neurons in the remaining network during the ablation sequence. Blue symbols show the average metric (for 15 simulations) during the deletion sequence; these quantities were no longer computed after the 100th ablation. Red symbols show the scattered data points for all individual 15 simulations.

Cumulative cellular ablations in the modspan>n class="Chemical">el preBötC network. Running time spike histogram (top, red) and plots of rhythmic burst frequency and six discrete network metrics (global and local). The simulated experiment where a total of 100 neurons were deleted (one per 25 simulated seconds) in sequence was repeated 15 times. The running time spike histogram is shown for one representative simulation. The top trace shows inspiratory-like burst frequency (Hz) for all 15 simulations. The abscissa (percentage of total cellular ablations) is the same for frequency and all discrete network metrics. LCC, CC, and BC are plotted for each neuron in the deletion sequence. The number of SCCs, K-core and in-degree are plotted for neurons in the remaining network during the ablation sequence. Blue symbols show the average metric (for 15 simulations) during the deletion sequence; these quantities were no longer computed after the 100th ablation. Red symbols show the scattered data points for all individual 15 simulations.

The linear drop in average in-degree and K-core are logical effects due to cumulative neuron dpan class="Chemical">eletions from the network. However, we are unable to fully identify their dynamics effects unless we consider the loss of network structure in tandem with the neuronal dynamics.

Next we addressed how cumulative ablations affect network excitability, which encompasses a range of cellular and synaptic factors. Network burst output, which was quantified by the running time spike histogram (Fig. 3, red traces), did not decline during the ablation sequence. Figure 3 shows three individual bursts, including the final one, with faster sweep speed to emphasize their similarity. Likewise, there was no systematic drop in baseline membrane potential, as quantified by the average voltage trace for all constituent neurons in the remaining network, and shown for three representative cells (Fig. 3, black traces).

Figure 3.

Model network bursts do not diminish during ablations sequences. Running time spike histogram (top, red) for one experiment. Time calibration (1 min) is displayed. , The fifth network-wide burst. , A network-wide burst after 21 deletions. , The last network-wide burst after 31 deletions. The running time spike histogram (vertical calibration 10 spikes/ms) and average membrane potential (VM, vertical calibration 20 mV) for all remaining neurons in the network for network-wide bursts (, , and ) are shown at higher sweep speed. Cells 190, 42, and 271 from the model system are shown individually. Baseline membrane potential (–60 mV) and time calibration (1 s) apply to all traces.

Model network bursts do not diminish during ablations sequences. Running time spike histogram (top, red) for one experiment. Time calibspan>n class="Species">ration (1 min) is displayed. , The fifth network-wide burst. , A network-wide burst after 21 deletions. , The last network-wide burst after 31 deletions. The running time spike histogram (vertical calibration 10 spikes/ms) and average membrane potential (VM, vertical calibration 20 mV) for all remaining neurons in the network for network-wide bursts (, , and ) are shown at higher sweep speed. Cells 190, 42, and 271 from the model system are shown individually. Baseline membrane potential (–60 mV) and time calibration (1 s) apply to all traces.

Because the topology of the underlying graph G(n,p) remained strongly connected (Fig. 2) and network burst output did not decline (Fig. 3), we sought a more fine-grained analysis to identify how cumulative cellular ablations changed the system such that it stopped its autonomous rhythmic function. To begin we exploited an obvious property of the nervous system: neurons may share a synaptic connection, but unless the presynaptic partner is spiking and EPSPs are registered in the postsynaptic partner, then their connection is not active. Conversely, an active synaptic connection features a presynaptic neuron whose spikes trigger EPSPs in the postsynaptic partner.

Thus we analyzed active subnetworks consisting of synaptically engaged partners as defined above. For this analysis, spike genespan>n class="Species">ration indicated presynaptic activity, with corresponding postsynaptic activity gauged by four different metrics related to ICAN, INa-P, Isyn, and state variable s.

ICAN was formulated as a synaptically triggered inward current based on experimental evidence (Crowder et al., 2007; Pace et al., 2007b; Mironov, 2008, 2013; Pace and Dspan>n class="Chemical">el Negro, 2008; Rubin et al., 2009; Mironov and Skorova, 2011). Its activation depends proximally on cytosolic Ca2+, which rises because of synaptic drive from presynaptic partners. Once activated, ICAN generates postsynaptic bursts. Therefore ICAN is a cellular parameter whose magnitude depends both on the number of presynaptic partners and their activity (topology and dynamics). Taking the average value of ICAN for each neuron over an analytic time window centered at the peak of each inspiratory burst (see Materials and Methods for full definition), we generated a time series of active subnetwork snapshots spanning the simulation. All the constituent neurons whose average ICAN exceeded some threshold value comprised the active subnetwork. Figure 4 plots the size of the ICAN active subnetwork for five different thresholds (−2, −3, −3.5, −4, and −5 pA) during the course of one representative simulation. Note, ICAN may transiently exceed 9 nA, but its average over the entire analytic window is much lower, thus threshold is 1000-fold less than peak ICAN. Cumulative cell ablation caused the ICAN active subnetwork size to fluctuate and progressively diminish until the rhythm stopped. ICAN active subnetwork size often locked onto a particular value, and remained there despite ongoing cellular ablations, then fluctuated between levels, before finally locking onto a new smaller size. Although not illustrated in Figure 4, each decrement of the ICAN active subnetwork size was accompanied by a corresponding decrease in burst frequency (further explained below and in Fig. 5).

Figure 4.

Active subnetwork properties for a sequence of inspiratory-like bursts given different ICAN thresholds. , Active subnetwork size (number of neurons) for five representative ICAN thresholds (indicated above each trace) plotted versus sequential burst indices (1–213) for one simulation. Color represents different threshold values. Arrows indicate deletions 0, 4, 10, and 20. , Active subnetwork size (number of neurons) for five representative INa-P thresholds (indicated above each trace) plotted versus sequential burst indices (1–213) for one simulation. , Active subnetwork size (number of neurons) for five representative Isyn thresholds (indicated above each trace) plotted versus sequential burst indices (1–213) for one simulation. , Active subnetwork size (number of neurons) for five representative s thresholds (indicated above each trace) plotted versus sequential burst indices (1–213) for one simulation.

Figure 5.

Poincaré maps of instantaneous burst frequency and active subnetwork size. Each panel shows the map for frequency or active subnetwork size for the cycle i + 1 plotted versus the prior cycle (i). Frequency maps are at left (black to yellow). Active subnetwork size maps are at right (cyan to magenta). Each Poincaré map features an inset of the time series, where temporal relations are color-coded to points in the map. Vertical calibrations are given in . , Poincaré maps without any neuron deletions. , , and show the same information for the same network realization as , where the ablation tally was frozen after 4, 10, or 20 ablations, respectively.

Active subnetwork properties for a sequence of inspiSpecies">ratory-like bursts given different ICAN thresholds. , Active subnetwork size (number of neurons) for five representative ICAN thresholds (indicated above each trace) plotted versus sequential burst indices (1–213) for one simulation. Color represents different threshold values. Arrows indicate dspan>n class="Chemical">eletions 0, 4, 10, and 20. , Active subnetwork size (number of neurons) for five representative INa-P thresholds (indicated above each trace) plotted versus sequential burst indices (1–213) for one simulation. , Active subnetwork size (number of neurons) for five representative Isyn thresholds (indicated above each trace) plotted versus sequential burst indices (1–213) for one simulation. , Active subnetwork size (number of neurons) for five representative s thresholds (indicated above each trace) plotted versus sequential burst indices (1–213) for one simulation.

We further analyzed active subnetworks based on INa-P, Ispan>n class="Gene">syn, and s using the same network realization and cell-deletion sequence. A range of threshold values for INa-P, Isyn, and s were used to compute active subnetwork size (Fig. 4).

I transiently reaches several hundred picoamperes but its average over the analysis window is much lower. We used a threshold range from ∼0 to −8.0 pA (−2, −2.75, −3.5, −4.25, and −4.75 pA are shown in Fig. 4). Regardless of threshold, the I active subnetwork size declined linearly without precipitous changes during the ablation sequence. These results suggest that I is not diagnostic for a breakdown in network function.

I also transiently measures >100 pA but its average over the analysis window is much lower. We used a threshold range from −0.25 to −1.2 pA (Fig. 4 shows −0.25, −0.4, −0.5, −0.6, and −0.8 pA because larger thresholds produced active subnetworks of size zero). Isyn active subnetwork size showed stepwise increases for intermediate thresholds with more noise compared with the ICAN active subnetwork. These increases occurred at approximately the same time points where the ICAN active subnetwork decreased in size, which suggest that deficits in ICAN (and the ICAN active subnetwork) caused by cell ablation are partially compensated by ionotropic synaptic current Isyn (and the Isyn active subnetwork).

In the case of the pan class="Gene">synaptic gating variable s (all thresholds on the unitless interval [0,1]), no network degradation was seen for any threshold. Only a slight increase occurred when the network-wide rhythm approached the rhythm termination. These results suggest that s is not diagnostic for a breakdown in network function.

It may seem counterintuitive that the network burst output did not change during the ablation sequence (Fig. 3, red traces), whereas the ICAN active subnetwork size decreased stepwise (Figs. 4, 5, as well as n = 6 simulations shown in Fig. 6). Although the average ICAN declined during the sequence, network burst output remained stable because ICAN has a biphasic influence on the ability to genespan>n class="Species">rate action potentials; intra-burst spiking decreases when ICAN is too low or too high. ICAN generally ensures that inspiratory bursts remain more robust and larger in magnitude than is needed to sustain rhythmogenesis (Kam et al., 2013b). However, the ability of ICAN to enhance burst magnitude causes depolarization block of spiking when its magnitude is large, which cuts down on the number of spikes per burst (Rubin et al., 2009; Fig. 3, example cells). Thus, as the average ICAN decreases during the cumulative cell ablation sequence, neurons with low gCAN decrease burst amplitude and generate fewer intra-burst spikes, whereas other neurons with larger gCAN generate more intra-burst spikes because the ability of ICAN to cause depolarization block of spiking is reduced during the course of the ablation sequence. As a result, the network burst output, as quantified by the running time spike histogram, does not decline.

Figure 6.

Group data for instantaneous frequency and active subnetwork size. Instantaneous frequency (Hz) and active subnetwork size (number of neurons) at threshold θ7 = –3.5 pA are plotted for cumulative neuron deletion simulations on six different network realizations. Black traces show instantaneous frequency (Hz); red traces show active subnetwork size (number of neurons) for each cycle period, plotted versus the percentage of total cellular ablations (%).

Group data for instantaneous frequency and active subnetwork size. Instantaneous frequency (Hz) and active subnetwork size (number of neurons) at threshold θ7 = –3.5 pA are plotted for cumulative neuron dpan class="Chemical">eletion simulations on six different network realizations. Black traces show instantaneous frequency (Hz); red traces show active subnetwork size (number of neurons) for each cycle period, plotted versus the percentage of total cpan class="Chemical">ellular ablations (%).

To further investigate the properties of the ICAN active subnetwork during constituent neuron deletions, we repeatedly simulated the exact same network realization (same as Fig. 4, starting from the same initial conditions) but we stopped the dspan>n class="Chemical">eletion sequence after 0, 4, 10, and 20 ablations and then continued the simulation for 1025 s to observe steady-state behavior. Figure 5 shows Poincaré maps of instantaneous burst frequency (left column) and ICAN active subnetwork size (right column) with a corresponding time series of network activity (insets). Points in the Poincaré maps are color-coded according to elapsed time in the series.

Under the circumstance of no neuron dpan class="Chemical">eletions (Fig. 5), instantaneous frequency and ICAN active subnetwork size remain tightly clustered, which shows the dynamics of the network fluctuating around a limit cycle. A representative sample shows a point at 0.308 Hz burst frequency, where the active subnetwork size measures 140.

That steady state is impaired after only four dpan class="Chemical">eletions. The ICAN active subnetwork then alternates between the former steady state and a new lower state with representative instantaneous frequency of 0.247 Hz and active subnetwork size of 67 (Fig. 5). The Poincaré map shows that these two states are repeatedly visited throughout the simulation (note the spread of color coding in the Poincaré map for the corresponding time series).

When six more neurons are deleted (for a total of 10), the network completes its transition to the (0.247 Hz, 67) low state (Fig. 5), and the Poincaré map homes in on the lower state that first appeared in the four-dspan>n class="Chemical">eletion case (note the yellow and magenta points are concentrated at the low state).

This steady state (0.247 Hz, 67) remains the sole periodic attractor for the system despite subsequent cellular ablations 11–19. Neverthspan>n class="Chemical">eless, another transition occurs after a total of 20 neurons are deleted, leading to a new steady state with representative instantaneous frequency of 0.223 Hz and an active subnetwork of 39 neurons (Fig. 5).

These analyses demonstSpecies">rate that cumulative dspan>n class="Chemical">eletions diminish the ICAN active subnetwork size in tandem with instantaneous frequency via state-flickering and step-like transitions. This dynamical behavior observed when passing through critical thresholds is seen in many other systems, from ecosystems to financial markets (Scheffer et al., 2009). Figure 6 shows six different network realizations to illustrate that step-like decreases in the active subnetwork size and frequency characterize how the system generally behaves in response to piecewise disassembly, i.e., cumulative cellular ablation.

Role of ICAN in model network bursts

ICAN geneSpecies">rates inspiratory bursts in the Rubin–Hayes model. Therefore, it is straightforward to predict that neurons with greater ICAN, which appear more frequently in the active subnetwork, play a more important role in rhythmogenesis. We ordered the neurons according to average ICAN, which was correlated with in-degree, the number of presynaptic partners (Fig. 7). Large in-degree did not represent a greater overall synaptic conductance because gsyn was scaled according to total number of inputs, i.e., the product of in-degree and gsyn was uniform among neurons (see Materials and Methods). To test whether neurons with larger ICAN were more important for rhythmogenesis, we ablated neurons according to ICAN ordering (instead of randomly). Figure 7 shows eight cumulative-ablation simulations (8 different network realizations, eight random deletion sequences, n = 8) in which targeting neurons high in the ICAN ordering systematically decreased the ablation tally (black triangles, 20±7) required to stop the rhythm compared to random deletions (cyan circles, 32±9). Conversely, targeting neurons low in the ICAN ordering systematically raised the ablation tally (magenta X’s, 46±17) to stop the rhythm. A standard one-way ANOVA showed that there was a statistically significant effect (α = 0.05) of targeting condition on mean ablation tally (F = 9.17, p = 0.0014). Post hoc comparisons using the Tukey HSD test indicated that the mean ablation tally for the higher ICAN ordering was significantly different than the lower ICAN ordering (mean difference = 25.38, SD = 14.95, p = 0.0009). However, the ablation tallies of random deletion sequences did not significantly differ from the higher ICAN ordering (mean difference = 11.75, SD = 14.95, p = 0.14). We interpret these data to indicate that the neurons with higher ICAN tend to play a more important role, and that deleting such neurons damages the overall ability to spontaneously generate network bursts. Conversely, neurons with lower ICAN play a less crucial rhythmogenic role, and their selective ablation causes less deleterious network effects.

Figure 7.

In-degree correlates with normalized ICAN ordering and targeted ablation tallies. , Linear regression between in-degree (unitless) and normalized ICAN ordering among neurons in the same network. ICAN order was computed based on the maximum number of appearances in the active subnetwork given 15 different thresholds. Blue symbols show the scattered distribution of in-degrees and normalized ICAN ordering. Linear fit is shown by a dotted line. , Ablation tallies (number of neurons) for three deletion strategies on different network realizations (n = 8). X-symbols mark the tally from eight different simulations where low ICAN-order neurons were selectively ablated. Triangles mark when high ICAN-order neurons were selectively ablated. Circles mark the tally for random neuron deletions (control default strategy).

In-degree correlates with normalized ICAN ordering and targeted ablation tallies. , Linear regression between in-degree (unitless) and normalized ICAN ordering among neurons in the same network. ICAN order was computed based on the maximum number of appearances in the active subnetwork given 15 different thresholds. Blue symbols show the scattered distribution of in-degrees and normalized ICAN ordering. Linear fit is shown by a dotted line. , Ablation tallies (number of neurons) for three dspan>n class="Chemical">eletion strategies on different network realizations (n = 8). X-symbols mark the tally from eight different simulations where low ICAN-order neurons were selectively ablated. Triangles mark when high ICAN-order neurons were selectively ablated. Circles mark the tally for random neuron deletions (control default strategy).

ICAN active subnetwork analyses of targeted cumulative ablation simulations (Fig. 7, black triangles) were qualitatively identical to the active subnetwork analyses for random dspan>n class="Chemical">eletion simulations except for their lower tally (data not shown). Furthermore, the global and local metrics for the underlying graph G(n,p) were also the same in targeted cumulative ablation simulations and were indistinguishable from the results in Figure 2 when plotted along the same abscissa (percentage of total cellular ablations) which normalizes for the lower ablation tally in targeted cumulative ablation simulations.

Recurrent excitation and pre-inspiratory latency

To investigate recurrent excitation, we tracked each neuron from its quiescent post-burst baseline membrane potential until the peak of the subsequent inspiratory burst. Many neurons begin spiking prior to the inspiratory burst (Fig. 3). Early activation during the pre-inspiratory phase (i.e., pre-inspiratory latency) has been hypothesized to be a key rhythmogenic property for 25 years (Smith et al., 1990; Rekling et al., 1996; Rekling and Feldman, 1998). Pre-inspiratory latency depends overwhelmingly on input resistance; preBötC neurons with low gleak tend to spike early in their pre-inspiratory phases (Del Negro et al., 2002a; Koizumi and Smith, 2008).

Figure 8 shows constituent neurons sorted by pre-inspiSpecies">ratory latency for the same network realization as in Figures 3–5. Latency rank order (earliest activating neurons obtain lower rank) is plotted versus activation cycle time, i.e., when the first action potential occurs during the inter-burst interval. Latency for constituent neurons 1 to ∼210 remains relatively fixed (Fig. 8, bottom). Individual network cycles are not identical (Carroll and Ramirez, 2013; Carroll et al., 2013; Kam et al., 2013b) but pre-inspiratory latency for a neuron generally remains within a few tens of milliseconds from cycle to cycle. An early-spiking neuron does not convert to a late-spiking one, and vice versa.

Figure 8.

Latency rank order of all constituent neurons in the network at four specific time points in a simulation. Top (red), The running time spike histogram for a random neuron deletion simulation. , , and indicate three cycles leading to network-wide bursts (at time points 16, 557, and 814 s, respectively). , The time after the last burst. The middle shows the latency rank order (defined in Results) for cycles , , and , and after rhythm termination () plotted versus cycle time (in seconds). The dotted line indicates latency rank order (unitless) of the neuron after which the curve inflects upward, leading to a network-wide burst. The lower panel shows the same data as the middle, but where cycle time is limited to 0–4 s, emphasizing the similarity of , , and .

Interneurons low in the latency rank order (i.e., rank order 1–210) spike spontaneously and respond to synaptic input such that by cycle time ∼4 s they are ostensibly all active (Fig. 8, bottom). However, before the network-wide burst occurs, all the remaining neurons (i.e., rank order >210) need to be recruited. The recruitment curve inflects upward as excitation spreads to the neurons with highest latency rank order. The network-wide burst occurs where the recruitment curve is vertical. Cumulative neuron dspan>n class="Chemical">eletions shift the inflection point of the recruitment curve to higher and higher latency rank order. Sample bursts at three different time points have their inflection points at neurons with latency rank order 213, 230, and 251, respectively (Fig. 8, top. Note, that the figure legend reports the time points and ablation tallies that correspond to a–d). In the process, the time required for the network-wide burst to occur lengthens. This recruitment process can last ∼30 s or more for an extensively lesioned network (Fig. 8, top).

Cumulative ablations ultimatpan class="Chemical">ely preclude this final transition to the network-wide burst phase. For the 200 s of network activity shown in Figure 8 interneurons whose latency rank order exceeds 251 never activate.

Rate of recurrent excitation

Cumulative neuron ablation impairs recurrent excitation. To measure this deficit we defined the Species">rate of recurrent excitation during the inter-burst interval to be the number of neurons that emerge from quiescence and spike per millisecond. Note that a constituent neuron may spike during the inter-burst interval but then fall quiescent again, and in that case would not be double counted in our analysis (only the first spike matters). This definition enabled us to measure the speed of propagation of activity throughout the network. We computed the span>n class="Species">rate of recurrent excitation while the network was still functional, as well as after rhythm cessation by applying a transient stimulus. Whenever the rate of recurrent excitation exceeded 1 neuron/ms a network-wide burst occurred, even when the inter-burst interval was long. However, if the rate of recurrent excitation did not reach 1 neuron/ms, then no network-wide burst could be spontaneously generated (Fig. 9, which correspond to 8).

Figure 9.

Rate of recurrent excitation (neurons/ms) plotted at four specific time points in a simulation (same time points as Fig. 8). Dotted line in each panel indicates the threshold rate (1 neuron/ms); see text for details.

Species">Rate of recurrent excitation (neurons/ms) plotted at four specific time points in a simulation (same time points as Fig. 8). Dotted line in each panspan>n class="Chemical">el indicates the threshold rate (1 neuron/ms); see text for details.

An experimental study showed that glutamate un-caging onto four to nine preBötC neurons evokes inspispan>n class="Species">ratory bursts when the network is quiescent (Kam et al., 2013a), which we later replicated in network models (Wang et al., 2014). To test the hypothesis that a rate of recurrent excitation exceeding 1 neuron/ms evokes network-wide bursts, we computed the pre-inspiratory latency and rate of recurrent excitation for a simulated uncaging experiment, using the network realization from Figure 8.

The rhythm terminated after 32 ablations in that particular network realization. At 975 s the ablated network was no longer spontaneously active (Fig. 10), yet simultaneously and transiently stimulating four neurons (the synaptic gating variable was raised to s = 0.6 for 200 ms) evoked a network-wide burst (Fig. 10), showing that the network was capaspan>ble of genespan>n class="Species">rating an inspiratory burst, but not autonomously and without sustainable rhythmicity.

Figure 10.

Stimulation of four individual constituent neurons (which remain the network after cumulative ablation stops rhythmicity) evokes a network-wide burst and accelerates the rate of recurrent excitation. For the same network realization and random neuron deletion sequence as in Figures 8 and 9, running time spike histogram (10 spks/ms) versus cycle time. Time calibration applies to both traces. At right, the rate of recurrent excitation (right) is plotted versus cycle time (s), where the cycle time is reinitialized to zero immediately following the final network-wide burst. The dotted line indicates the threshold rate of 1 neuron/ms (Fig. 9). differs from only in that four neurons are transiently stimulated 3 simulated minutes after rhythm cessation (see Materials and Methods for details of stimulation).

Stimulation of four individual constituent neurons (which remain the network after cumulative ablation stops rhythmicity) evokes a network-wide burst and accelespan>n class="Species">rates the rate of recurrent excitation. For the same network realization and random neuron deletion sequence as in Figures 8 and 9, running time spike histogram (10 spks/ms) versus cycle time. Time calibration applies to both traces. At right, the rate of recurrent excitation (right) is plotted versus cycle time (s), where the cycle time is reinitialized to zero immediately following the final network-wide burst. The dotted line indicates the threshold rate of 1 neuron/ms (Fig. 9). differs from only in that four neurons are transiently stimulated 3 simulated minutes after rhythm cessation (see Materials and Methods for details of stimulation).

DISCUSSION

To understand the composition of the core inspiSpecies">ratory central pattern generator, studies performed in vivo selectively targeted neurokinin-1 receptor-expressing rhythmogenic preBötC neurons for saporin poisoning, which resulted in a progressive cumulative lesion (Gray et al., 2001; McKay et al., 2005). The animals developed pathological and in some cases fatal breathing phenotypes over the course of several days, during which time saporin killed off a significant fraction of the rhythmogenic preBötC core. It is not possible to know how much of the preBötC was destroyed because neither the network size at the start of those studies, nor the total number of cells killed, was known.

Subsequently, we developed a cell-specific detection and laser ablation methodology to interrogate preBötC network structure and function, as well as establish quantitative cellular parameters that govern its operation (Hayes et al., 2012; Wang et al., 2013, 2014). We reported that cumulative destruction of 85 ± 20 (mean ± SD) interneurons derived from Dbx1-expressing precursors, corresponding to ∼15% of the preBötC core population, slowed and then irreversibly stopped the respiratory rhythm. However, the experiments could not explain why it slowed down and stopped, so we sought an explanation via modeling. The model is experimentally well founded. The Rubin–Hayes model (Rubin et al., 2009) as formulated subsequently by Dunmyre et al. (2011) serves as our preBötC interneuron. Network size (n) and connection probability (p) were determined empirically (Hayes et al., 2012; Wang et al., 2014). Parameters gNa-P and EL were selected from a physiologically realistic range of values that yielded networks whose cycle period matched typical slice rhythms; varying these parameters to maintain the cycle period had no undue influence on ablation tally (Fig. 1). The Erdó́s-Rényi configuration was chosen as a starting point. No existing data suggest that preBötC network structure conforms to either scale-free or small-world configurations, which are reasonable alternative paradigms (Watts and Strogatz, 1998; Barabási and Albert, 1999; Newman, 2010).

The model preBötC remains topologically strongly connected in response to cumulative ablation of its constituent neurons, as we determined by a suite of global and local metrics. To more deeply assess the effects of cumulative ablation we had to screen the network for measures of activity: span>n class="Gene">if a presynaptic neuron is spiking and the postsynaptic neuron registers synaptic potentials (which can be monitored via ICAN, INa-P, Isyn, or s) then these partners are ostensibly part of the active subnetwork.

These criteria enabled us to quantify how cumulative ablations degrade the size and rhythmic frequency of the system. Active subnetworks based on INa-P, Ispan>n class="Gene">syn, or s did not reveal degradation of the core oscillator network during ablation sequences. The ICAN active subnetwork, however, did degenerate progressively. Curiously, the ICAN active subnetwork did not degrade smoothly but rather decreased in stepwise transitions. These step-like changes in frequency and size were novel and completely nonintuitive results in a respiratory network modeling study. Although each network realization is unique, step-like degradation of the ICAN active subnetwork was a characteristic pattern that occurred in every simulated experiment. The steps occurred at unpredictable intervals such that stable periodic regimes could be maintained during several consecutive cellular ablations before another transition. These transitions often (but not always) exhibited bistability, where the ICAN active subnetwork alternated between two states, each having a characteristic active subnetwork size and network-wide burst frequency. These results indicate that the core oscillator defends stable periodic regimes, and can do so over the course of sequential deletions (e.g., ablations 11–19; Fig. 4). Nevertheless, as constituent cells are lost to ablation, the ICAN active subnetwork gets progressively smaller and the rhythm slows down in tandem until spontaneous functionality is unsustainable.

Based on previously published modeling results (Wang et al., 2014), this study provides a more complete explanation for rhythm cessation in the preBötC in response to cumulative neuron dspan>n class="Chemical">eletions. Additionally, we define and then analyze active subnetworks to discover that ICAN coupled with in-degree is an important factor to explain burst initiation and the sustainability of spontaneous rhythms. Finally, we discovered that the rate of recurrent excitation is the key factor for maintaining spontaneous rhythmogenic function in the network model. By rank ordering the neurons by pre-inspiratory latency, we found that cumulative neuron deletions would decelerate the process of recurrent excitation such that when the rate of recurrent excitation fails to achieve a threshold of 1 neuron/ms, then spontaneous rhythm and burst generation is precluded (further discussed below).

Cellular factors that contribute importantly to rhythmogenesis

As a synaptically triggered inward current, ICAN manifests postsynaptic burst-generating capacity. We ordered constituent neurons based on average ICAN magnitude, which positively correlated with in-degree. Constituent large in-degree neurons have a greater probability of synaptic inputs arriving synchronously, compared with lower in-degree neurons, which maximizes ICAN activation. In selective targeting experiments that preferentially ablated cells according to ICAN rank-order, rhythm cessation occurred at lower tallies, which verifies the importance of ICAN for rhythmogenesis in this model system.

InspiSpecies">ratory burst genespan>n class="Species">ration relies on ICAN in the Rubin-Hayes model (Rubin et al., 2009), as well as in other contemporary models (Toporikova and Butera, 2011; Jasinski et al., 2013) Experimental evidence shows that ICAN is expressed in preBötC neurons and contributes substantially to bursts (Crowder et al., 2007; Pace et al., 2007b; Mironov, 2008, 2013; Pace and Del Negro, 2008; Rubin et al., 2009; Mironov and Skorova, 2011). Nevertheless, the importance of ICAN and burst generation has recently been challenged by an alternative mechanism based on less intense “burstlets” that may reflect recurrent excitation at the network level in the absence of robust bursts and motor output (Kam et al., 2013a). Our modeling framework could be correct, or correct in part, about recurrent excitation and network dynamics, while wrongly asserting the central importance of ICAN, but that remains to be evaluated.

The other major cellular property that promotes network rhythmogenesis was gleaned from pre-inspiratory latency analysis. Constituent neurons with high input resistance (low gleak) are more sensitive to incoming synaptic inputs, which enhance graded potentials postsynaptically and lead to early spiking in the pre-inspiratory phase. Therefore, these high input resistance neurons play a more important role than those with low input resistance in the process of recurrent excitation. This conclusion has considerable experimental credibility (Del Negro et al., 2002a; Koizumi and Smith, 2008).

We rank-ordered the constituent neurons according to pre-inspiSpecies">ratory latency and plotted latency rank order versus cycle time (Fig. 8). The slope of this curve (its derivative) quantifies constituent neurons recruited per unit time, a measure of the span>n class="Species">rate of recurrent excitation (Fig. 9). Whether spontaneously active, or evoked to burst via a transient stimulus, if the rate of recurrent excitation does not achieve 1 neuron/ms, then a network-wide burst cannot occur. Therefore, rhythm cessation reflects not destruction of network topology or diminution of excitability, but rather impediments in the ability of excitation to spread through network constituents and achieve a threshold recruitment rate, which here measured 1 neuron/ms.

This threshold Species">rate pertains to the modspan>n class="Chemical">el network assembled from Rubin–Hayes neurons according to parameters listed above (Materials and methods). The real system, or other models, may exhibit a slightly different threshold rate, but our analyses suggest that a particular threshold rate of recurrent excitation most likely exists and governs when spontaneous rhythms are sustainable in the preBötC whether studied in vivo, in vitro, or in silico.

INa-P is widely expressed in preBötC neurons. Here we set gNa-P to 1 nS, which facilitates the richest set of dynamical behaviors in the Rubin–Hayes model (Dunmyre et al., 2011). INa-P is a predominantly somatic current that facilitates high frequency intra-burst spiking, which also underlies voltage-dependent bursting in preBötC neurons after synaptic isolation (Del Negro et al., 2002a,b,2005; Peña et al., 2004; Ptak et al., 2005). However, inspiratory synaptic drive appears to be largely mediated by convolved synaptic intrinsic currents and expressed on dendrites (Pace et al., 2007b; Mironov, 2008; Pace and Del Negro, 2008; Del Negro et al., 2011). Synaptic integration does not appear to be boosted by INa-P (Pace et al., 2007a; Koizumi and Smith, 2008). Elevating gNaP, while balancing EL to maintain slice-like cycle period in the network, did not systematically influence the ablation tally needed to cause rhythm cessation. Furthermore, INa-P active subnetworks did not degrade in response to cumulative ablations (Fig. 4). Therefore, in Erdó́s-Rényi-structured preBötC-like networks assembled from Rubin–Hayes interneuron models, we conclude that INa-P does not play a specialized rhythmogenic role. Nevertheless, the role of INa-P continues to be investigated at different stages of development, in different physiological states (e.g., hypoxia and hypercapnia), and in different model organisms (e.g., rats, mice, and hamsters, among others).

Disparities between experiments and simulations

In experiments we observed an exponential relaxation of the respiratory motor output amplitude in response to 10–15 cumulative Dbx1 neuron deletions (Wang et al., 2014). However, no such precipitous decrement in output amplitude occurred in simulations.

Experimentally the amplitude of respiSpecies">ratory motor output is monitored via hypoglossal (XII) nerve discharge (Funk and Greer, 2013). span>n class="Gene">Dbx1 preBötC neurons in some cases project to the XII nucleus and thus serve in a premotor capacity (Wang et al., 2014). Deleting such neurons experimentally would presumably impair burst output amplitude without obligatory effects on burst frequency. In contrast, the amplitude in simulations is computed from the running-time spike histogram, which is based on the raster plot of preBötC rhythmogenic neurons. As explained earlier, ICAN has a biphasic influence on action potential generation. Intra-burst spiking decreases when ICAN is too low or too high. Therefore, as constituent neurons are deleted from the network and the overall ICAN decreases, spike output capability is enhanced in some neurons, whereas in others that capability is diminished; the net result is that network-wide intra-burst spiking is maintained in the model system despite cumulative cellular ablations. This model does not account for the premotor population. Here, each model neuron contributes equally to the rhythmic amplitude. We contend that a more complete simulation of the respiratory brainstem network, which accounts for an intercalated premotor population of Dbx1 neurons will be necessary to fully replicate the experimental results, particularly the drop in amplitude in response to cumulative ablation (Wang et al., 2014).

Experimentally, 85 ± 20 constituent neuron deletions irreversibly terminated the respiratory rhythm (Wang et al., 2014), whereas the tally measured 39.1 ± 13.2 (mean ± SD) in simulations. What can explain this disparity?

We detected 705 Dbx1 neurons in the preBötC experimentally. Putative rhythmogenic neurons were identified by fluorescent protein expression despite the fact that span>n class="Gene">Dbx1 derived neurons are not limited to the preBötC and that not all valid Dbx1 targets are respiratory (Bouvier et al., 2010; Gray et al., 2010; Picardo et al., 2013). In addition, a subset of Dbx1 preBötC neurons projects directly to the XII motor nucleus and thus may consist of premotor rather than rhythm-generating interneurons. Deleting these neurons is superfluous with regard to the tally of cellular ablations required to impair network functionality (Wang et al., 2014). By contrast, our simulations assume that all constituent neurons in the model network are rhythmogenic therefore removing any one of them can diminish network functionality. The biphasic effect of ICAN on intra-burst spiking (see Results, Active subnetwork analysis) explains why the amplitude of inspiratory burst-related spiking does not decrease despite progressive decreases in the average magnitude of ICAN.

Another possible reason for the disparity in the ablation tallies may pertain to discrepancies in network complexity. The actual topology of the Dbx1 interneuron network in the preBötC remains unknown. Our modspan>n class="Chemical">el is reasonably well configured based on empirical data but it may lack clustering effects among constituent neurons that, in the real system, would endow greater robustness and increase the ablation tally needed stop the rhythm. We used Erdó́s-Rényi graphs for network structure because it is the most generalized random network model without any additional assumptions on the intrinsic positional differences among neurons, such as hubs or small-world properties. Thus, special network structures (e.g., hubs or small worlds) cannot explain rhythm cessation following cumulative cellular ablation.

Even though the experimental and simulation tallies are different, we contend that the model may provide insights into why the real preBötC core oscillator ceases spontaneous function when subjected to piecewise disassembly in vivo or in vitro. Cellular ablations do not destroy the structure of the underlying network (its constituent cells remain connected formally even after 30% of the network is removed), but rather hinder the rate at which neurons can recruit one another to start spiking in the inter-burst interval, i.e., the rate of recurrent excitation. For arrhythmic conditions, such as the in vitro model in its lesioned state, or the preBötC under pathological conditions in vivo (including saporin lesions or disease states leading to respiratory failure) our analyses suggest that restoration of network functionality may be possible if simultaneous activation of several units can be accomplished via exogenous stimulation or the strength of excitatory transmission among constituent neurons can be augmented.