Mechanisms Leading to Rhythm Cessation in the Respiratory PreBötzinger Complex Due to Piecewise Cumulative Neuronal Deletions1,2,3

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


Introduction
Inspiratory breathing movements originate within the preBötzinger complex (preBötC) of the lower medulla Feldman and Del Negro, 2006;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., 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 pre-Bö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.
Inspiratory neural bursts in the preBötC depend on excitatory synaptic transmission 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-offunction 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.

Synaptic inputs and calcium dynamics
Synaptically triggered increases in cytosolic Ca 2ϩ directly activate I CAN . We coupled the synaptic variable s to the Ca 2ϩ equation. The synaptic variable s represents both ionotropic and metabotropic glutamatergic receptor (mGluR) activation. Cytosolic Ca 2ϩ changes are attributable to influx through voltage-gated Ca 2ϩ channels evoked by AMPA receptor-mediated depolarization as well as group I mGluRs that evoke intracellular Ca 2ϩ 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 Ca 2ϩ , 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 Ca 2ϩ in all the postsynaptic neurons to which it projects. The parameter k IP3 governs how much Ca 2ϩ increases per unit increment in s. Variable s also appears in the chord-conductance equation for I syn , where it controls AMPA receptor gating.

Parameters
The model maintains rhythmic function over substantial variation in g CAN , g NaϪP , k IP3 , g syn , r pump , and k Ca from nominal baseline values when I pump is present. Initial values for the membrane potential and gating variables were set constant. The parameters g CAN , g leak , and overall g syn were randomly assigned from Gaussian distributions given the average and SD listed above. Then, g syn 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 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 generated 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 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 capacity, and peak performance of 21.2 teraflops. We used a Runge-Kutta fourth-order numerical integration 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 respiratory-related hypoglossal nerve (XII) root discharge in vitro 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 model 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 ij ϭ 1 if and only if there is a directed link (synaptic connection) from neuron i to neuron j; otherwise, A ij ϭ 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 i with k i links, the local cluster coefficient is defined as follows: where N i is the neighborhood of v i , (the subgraph formed by all the nodes v i connects to, that is, all the outneighbors of v i ; Watts and Strogatz, 1998). The numerator is the number of actual connections within N i , whereas the denominator is the number of connections if N i is a complete directed graph.

Closeness centrality
The farness of node v i is defined as the sum of shortest paths from v i to every other node reachable by v i along directed paths in the network. Closeness of v i is the reciprocal of the 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), st ͑v͒ is the number of shortest paths between s and t through v, while st 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. Nonetheless, 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 synaptic connection and an outgoing 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 d i out ϭ ͚

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). If a structural connection from neuron i to j (i.e., A ij ϭ 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 I CAN , which influences burst generation in the Rubin-Hayes model. The model also includes inward currents I Na-P and I syn . We examined active subnetworks based on I CAN , I Na-P , I syn , and s (which gates I syn , and indirectly, I CAN ).
Each simulation consists of N bursts before rhythm cessation, with peak times (t x ) of the inspiratory-like bursts (t 1 , t 2 , . . ., t N ) and the periods (T x ) between them (T 1 , T 2 , . . ., T N-1 ). We define an analytic window equal to min͑T xϪ1 , T x ͒, where xʦ͓2, N Ϫ 1͔, centered at t x . The last time window is equal to T N-1 , centered at t N . For each neuron at each time step, we compute the average I CAN (in units of pA) over its corresponding analytic window. Although peak I CAN transiently reaches 9 -15 nA, the average I CAN over the duration of the analytic window is much smaller, ranging from 0 to 10 pA. Then, we identify subsets of neurons whose I CAN 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 I CAN 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). E L and g Na-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 I Na-P (Dunmyre et al., 2011). At the start of each simulation (given n, p, E L , and g Na-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 model lacks, could explain, at least in part, this modelexperiment disparity (Wang et al., 2014), we did not previously investigate whether the excitability parameter E L or the conductance g Na-P influenced the ablation tally. To do so here, we first tested how E L and g Na-P influence network behavior. Using the same network realization, i.e., the same underlying G(n,p) structure, we simulated networks with E L spanning from Ϫ60.6 to Ϫ62.5 mV (E L outside this range is not physiologically realistic) and g Na-P spanning from 1 to 1.5 nS (g Na-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 E L or g Na-P slowed down the rhythm, such that for some (E L , g Na-P ) pairs the rhythm stopped ( Fig. 1, black squares) and for some the period would exceed 10 s ( Fig. 1, grey squares), whereas elevating E L or g Na-P had the opposite effect (it speeds up the rhythm) such that for some (E L , g Na-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 (E L , g Na-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 (E L , g Na-P ) combinations that yield networks with cycle period in the range 3.5-5 s are equally sensitive to cumulative cellular ablation.
Given that various (E L , g Na-P ) pairs produce realistic network rhythms, and the trade-off between E L and g Na-P results in commensurate ablation tallies, we used g Na-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 E L ϭ Ϫ61.46 mV, which straddles the two blocks on the diagonal line in Figure 1.

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.
The linear drop in average in-degree and K-core are logical effects due to cumulative neuron deletions 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.

Active subnetwork analysis
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).
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 finegrained 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 generation indicated presynaptic activity, with corresponding postsynaptic activity gauged by four different metrics related to I CAN , I Na-P , I syn , and state variable s. I CAN was formulated as a synaptically triggered inward current based on experimental evidence (Crowder et al., 2007;Pace et al., 2007b;Mironov, 2008Mironov, , 2013Pace and Del Negro, 2008;Rubin et al., 2009;Mironov and Skorova, 2011). Its activation depends proximally on cytosolic Ca 2ϩ , which rises because of synaptic drive from presynaptic partners. Once activated, I CAN generates postsynaptic bursts. Therefore I CAN is a cellular parameter whose magnitude depends both on the number of presynaptic partners and their activity (topology and dynamics). Tak- ing the average value of I CAN 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 I CAN exceeded some threshold value comprised the active subnetwork. Figure 4A plots the size of the I CAN active subnetwork for five different thresholds (Ϫ2, Ϫ3, Ϫ3.5, Ϫ4, and Ϫ5 pA) during the course of one representative simulation. Note, I CAN may transiently exceed 9 nA, but its average over the entire analytic window is much lower, thus threshold is 1000-fold less than peak I CAN . Cumulative cell ablation caused the I CAN active subnetwork size to fluctuate and progressively diminish until the rhythm stopped. I CAN 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 I CAN active subnetwork size was accompanied by a corresponding decrease in burst frequency (further explained below and in Fig. 5).
We further analyzed active subnetworks based on I Na-P , I syn , and s using the same network realization and celldeletion sequence. A range of threshold values for I Na-P , I syn , and s were used to compute active subnetwork size (Fig. 4B-D).
I Na-P 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. 4B). Regardless of threshold, the I Na-P active subnetwork size declined linearly without precipitous changes during the ablation sequence. These results suggest that I Na-P is not diagnostic for a breakdown in network function.
I syn 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. 4C shows Ϫ0.25, Ϫ0.4, Ϫ0.5, Ϫ0.6, and Ϫ0.8 pA because larger thresholds produced active subnetworks of size zero). I syn active subnetwork size showed stepwise increases for intermediate thresholds with more noise compared with the I CAN active subnetwork. These increases occurred at approximately the same time points where the I CAN active subnetwork decreased in size, which suggest that deficits in I CAN (and the I CAN active subnetwork) caused by cell ablation are partially compensated by ionotropic synaptic current I syn (and the I syn active subnetwork).
In the case of the 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 I CAN active subnetwork size decreased stepwise (Figs. 4A, 5, as well as n ϭ 6 simulations shown in Fig. 6). Although the average I CAN declined during the sequence, network burst output remained stable because I CAN has a biphasic influence on the ability to generate action potentials; intra-burst spiking decreases when I CAN is too low or too high. I CAN 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 I CAN 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 I CAN New Research decreases during the cumulative cell ablation sequence, neurons with low g CAN decrease burst amplitude and generate fewer intra-burst spikes, whereas other neurons with larger g CAN generate more intra-burst spikes because the ability of I CAN 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.
To further investigate the properties of the I CAN 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 deletion 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 I CAN 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 deletions (Fig.  5A), instantaneous frequency and I CAN 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 deletions. The I CAN 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. 5B). 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. 5C), and the Poincaré map homes in on the lower state that first appeared in the four-deletion 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. Nevertheless, 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.  5D).
These analyses demonstrate that cumulative deletions diminish the I CAN 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 I CAN in model network bursts
I CAN generates inspiratory bursts in the Rubin-Hayes model. Therefore, it is straightforward to predict that neurons with greater I CAN , which appear more frequently in the active subnetwork, play a more important role in rhythmogenesis. We ordered the neurons according to average I CAN , which was correlated with in-degree, the number of presynaptic partners (Fig. 7A). Large in-degree did not represent a greater overall synaptic conductance because g syn was scaled according to total number of inputs, i.e., the product of in-degree and g syn was uniform among neurons (see Materials and Methods). To test whether neurons with larger I CAN were more important for rhythmogenesis, we ablated neurons according to I CAN ordering (instead of randomly). Figure 7B shows eight cumulative-ablation simulations (8 different network realizations, eight random deletion sequences, n ϭ 8) in which targeting neurons high in the I CAN 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 I CAN 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 I CAN ordering was significantly different than the lower I CAN 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 I CAN ordering (mean difference ϭ 11.75, SD ϭ 14.95, p ϭ 0.14). We interpret these data to indicate that the neurons with higher I CAN 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 I CAN play a less crucial rhythmogenic role, and their selective ablation causes less deleterious network effects.
I CAN active subnetwork analyses of targeted cumulative ablation simulations (Fig. 7B, black triangles) were qualitatively identical to the active subnetwork analyses for random deletion 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 g leak 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 preinspiratory 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. 8a,b-d, 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.
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. 8a,b-d, 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 deletions 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. 8a-c, 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 re- cruitment process can last ϳ30 s or more for an extensively lesioned network (Fig. 8c, top).
Cumulative ablations ultimately preclude this final transition to the network-wide burst phase. For the 200 s of network activity shown in Figure 8d 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 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 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. 9a-d, which correspond to 8a-d).
An experimental study showed that glutamate uncaging onto four to nine preBötC neurons evokes inspiratory 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. 10A), 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. 10B), showing that the network was capable of generating an inspiratory burst, but not autonomously and without sustainable rhythmicity.

DISCUSSION
To understand the composition of the core inspiratory 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 pre-Bö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., , 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 g Na-P and E L 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: if a presynaptic neuron is spiking and the postsynaptic neuron registers synaptic potentials (which can be monitored via I CAN , I Na-P , I syn , 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 I Na-P , I syn , or s did not reveal degradation of the core oscillator network during ablation sequences. The I CAN active subnetwork, however, did degenerate progressively. Curiously, the I CAN 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 I CAN 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 I CAN 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. 4A). Nevertheless, as constituent cells are lost to ablation, the I CAN 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 deletions. Additionally, we define and then analyze active subnetworks to discover that I CAN 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, I CAN manifests postsynaptic burst-generating capacity. We ordered constituent neurons based on average I CAN magnitude, which positively correlated with in-degree. Constituent large in-degree neurons have a greater probability of synaptic inputs arriving synchronously, compared with lower indegree neurons, which maximizes I CAN activation. In selective targeting experiments that preferentially ablated cells according to I CAN rank-order, rhythm cessation occurred at lower tallies, which verifies the importance of I CAN for rhythmogenesis in this model system.
Inspiratory burst generation relies on I CAN 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 I CAN is expressed in preBötC neurons and contributes substantially to bursts (Crowder et al., 2007;Pace et al., 2007b;Mironov, 2008Mironov, , 2013Pace and Del Negro, 2008;Rubin et al., 2009;Mironov and Skorova, 2011). Nevertheless, the importance of I CAN 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 I CAN , 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 g leak ) 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-inspiratory 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 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 rate pertains to the model 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. I Na-P is widely expressed in preBötC neurons. Here we set g Na-P to 1 nS, which facilitates the richest set of dynamical behaviors in the Rubin-Hayes model (Dunmyre et al., 2011). I Na-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,2005Peñ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 I Na-P (Pace et al., 2007a;Koizumi and Smith, 2008). Elevating g NaP, while balancing E L to maintain slice-like cycle period in the network, did not systematically influence the ablation tally needed to cause rhythm cessation. Furthermore, I Na-P active subnetworks did not degrade in response to cumulative ablations (Fig. 4B). Therefore, in Erdó´s-Rényistructured preBötC-like networks assembled from Rubin-Hayes interneuron models, we conclude that I Na-P does not play a specialized rhythmogenic role. Nevertheless, the role of I Na-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). How-  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 9A,B, 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). B differs from A only in that four neurons are transiently stimulated 3 simulated minutes after rhythm cessation (see Materials and Methods for details of stimulation). ever, no such precipitous decrement in output amplitude occurred in simulations.
Experimentally the amplitude of respiratory motor output is monitored via hypoglossal (XII) nerve discharge (Funk and Greer, 2013). 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, I CAN has a biphasic influence on action potential generation. Intraburst spiking decreases when I CAN is too low or too high. Therefore, as constituent neurons are deleted from the network and the overall I CAN 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).
We detected 705 Dbx1 neurons in the preBötC experimentally. Putative rhythmogenic neurons were identified by fluorescent protein expression despite the fact that 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 I CAN on intra-burst spiking (see Results, Active subnetwork analysis) explains why the amplitude of inspiratory burstrelated spiking does not decrease despite progressive decreases in the average magnitude of I CAN .
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 model is reasonably well configured based on empirical data but it may lack clustering effects among constit-uent 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.