Effects of persistent sodium current blockade in respiratory circuits depend on the pharmacological mechanism of action and network dynamics

The mechanism(s) of action of most commonly used pharmacological blockers of voltage-gated ion channels are well understood; however, this knowledge is rarely considered when interpreting experimental data. Effects of blockade are often assumed to be equivalent, regardless of the mechanism of the blocker involved. Using computer simulations, we demonstrate that this assumption may not always be correct. We simulate the blockade of a persistent sodium current (INaP), proposed to underlie rhythm generation in pre-Bötzinger complex (pre-BötC) respiratory neurons, via two distinct pharmacological mechanisms: (1) pore obstruction mediated by tetrodotoxin and (2) altered inactivation dynamics mediated by riluzole. The reported effects of experimental application of tetrodotoxin and riluzole in respiratory circuits are diverse and seemingly contradictory and have led to considerable debate within the field as to the specific role of INaP in respiratory circuits. The results of our simulations match a wide array of experimental data spanning from the level of isolated pre-BötC neurons to the level of the intact respiratory network and also generate a series of experimentally testable predictions. Specifically, in this study we: (1) provide a mechanistic explanation for seemingly contradictory experimental results from in vitro studies of INaP block, (2) show that the effects of INaP block in in vitro preparations are not necessarily equivalent to those in more intact preparations, (3) demonstrate and explain why riluzole application may fail to effectively block INaP in the intact respiratory network, and (4) derive the prediction that effective block of INaP by low concentration tetrodotoxin will stop respiratory rhythm generation in the intact respiratory network. These simulations support a critical role for INaP in respiratory rhythmogenesis in vivo and illustrate the importance of considering mechanism when interpreting and simulating data relating to pharmacological blockade.


Introduction
Pharmacological compounds that selectively block voltage-gated ion channels are a fundamental tool in neuroscience. Much of our current theoretical understanding of the roles of various ion channels derives from the interpretation of data from experiments dependent on pharmacological manipulations. The mechanisms of action for many of the most commonly used pharmaceutical blockers of voltage-gated ion channels are relatively well understood and fall into one of three mechanistic categories: (1) pore obstruction, (2) shift in activation/inactivation curves, or, less commonly, (3) alteration of ion selectivity [1,2]. Despite this knowledge, the specific mechanisms of blockade are rarely considered when interpreting or simulating experimental data. Generally it is assumed that selective blockade of an ion channel has the same functional implication regardless of the mechanism involved. In this theoretical study, we demonstrate ways that this assumption can break down, with different blockade mechanisms differentially impacting neuronal and circuit activity.
To illustrate this idea we simulated blockade of a persistent sodium current (I NaP ) in the respiratory pre-Bötzinger complex (pre-BötC) via two commonly used sodium channel blockers with distinct mechanisms of action: tetrodotoxin (TTX) and riluzole (RZ). TTX directly obstructs the Na+ pore [3], whereas RZ shifts I NaP inactivation in the hyperpolarizing direction [4][5][6]. At low concentrations, both TTX (� 20 nM) and RZ (� 20 μM) have been shown to selectively block I NaP over the fast action-potential generating sodium current (I Na ) [7][8][9][10]. A caveat to RZ blockade, however, is that RZ has been shown to inhibit excitatory synaptic transmission at concentrations that affect I NaP [10,11].
In respiratory circuits, in vitro blocking studies using TTX or RZ have suggested that I NaP is critical for intrinsic bursting in pacemaker neurons and network rhythm generation in the isolated pre-BötC [9,12,13]. These results have led to the hypothesis that I NaP may be a necessary component of rhythmogenesis in respiratory circuits [14][15][16]. This hypothesis, however, has fallen out of favor due to the observation that I NaP block by RZ fails to stop respiratory rhythms in intact preparations [13]. This conclusion is dependent on the assumption that RZ effectively blocks I NaP in an in vivo setting.
This study predicts that after blockade via RZ, I NaP can be reactivated by transient hyperpolarizing perturbations due to the specific mechanism of action of RZ. Simulated TTX blockade of I NaP does not yield this same effect. In intact preparations, the pre-BötC receives strong inhibition during the interburst interval [13,[17][18][19][20], which may allow I NaP to recover from inactivation even after RZ application. Consistent with this idea, our simulations of the intact respiratory network predict that RZ will fail to effectively block I NaP , while complete block of I NaP by experimental application of low concentration TTX (� 20 nM) within the pre-BötC will abolish the respiratory rhythm in the intact respiratory network. Therefore, the failure of RZ to stop respiratory rhythms in intact experimental preparations is not sufficient to rule out a central role for I NaP in respiratory rhythm generation, which our study supports.
More generally, these simulations illustrate the importance of considering the mechanism of action when interpreting and simulating experimental data from pharmacological blocking studies.

Results
To illustrate the difference in effects that can arise through blockade of the same current with pharmacological agents acting through different biophysical mechanisms, we focused on the blockade of I NaP via RZ and TTX in respiratory neurons of the pre-BötC. Experiments have clearly established the presence of I NaP in neurons of this type and have suggested that neurons exhibiting I NaP -dependent intrinsic bursting capabilities play a critical role in rhythm generation in the pre-BötC [9,12,14]. Therefore, we first reconstructed a neuron model capable of generating I NaP -dependent intrinsic bursting [14]. For this set of simulations, the model represents a single pre-BötC neuron that only receives tonic excitatory synaptic drive with intensity set by the constant conductance parameter g Tonic . Depending on their level of excitability (controlled by g Tonic ), these neurons can exhibit three distinct patterns of activity: silent, bursting and tonic spiking ( Fig 1A). In an intrinsic bursting regime in this model, I NaP is required to depolarize neurons above the spiking threshold. During the spiking within a burst, however, I NaP begins to inactivate, which leads to hyperpolarization and burst termination. In the subsequent silent period or interburst interval, the membrane potential V m is hyperpolarized, which allows for the voltage-dependent recovery of I NaP from inactivation. After sufficient recovery, sub-threshold activation of I NaP again depolarizes V m above the spiking threshold, causing burst initiation and completion of one burst cycle (Fig 1B).

Simulated TTX and RZ blockade of I NaP in steady-state conditions
Bursting in this model is dependent on I NaP , therefore, we first characterized the effects of simulated TTX and RZ blockade on I NaP under steady-state conditions (Fig 2). Since TTX directly obstructs the pores of sodium channels, TTX blockade was simulated by reducing the conductance g NaP of I NaP . Since TTX only affects conductance, the steady-state activation and inactivation parameters of I NaP were not varied ( Fig 2A1). We observed that the steady-state current is reduced proportionally to the decrease of g NaP (Fig 2A2). In contrast, RZ shifts voltagedependent inactivation in the hyperpolarizing direction. Therefore, RZ blockade was simulated by decreasing the half-inactivation parameter h 1/2 ( Fig 2B1). We found that, similar to TTX, simulated RZ application effectively blocks I NaP in steady-state conditions ( Fig 2B2). Interestingly, although TTX and RZ blockade of I NaP work through distinct mechanisms, the effects on the peak steady-state current are remarkably similar ( Fig 2C).

Simulated TTX and RZ blockade of I NaP abolish intrinsic bursting in an isolated model pre-BötC neuron
The activity pattern exhibited by a neuron capable of I NaP -dependent bursting is determined by its excitability. Driving the neuron from low to high excitability will drive the transition from silent to bursting and from bursting to tonic spiking regimes of activity. Therefore, we investigated the effects of simulated TTX and RZ blockade of I NaP , represented by reduction of g NaP or h 1/2 , respectively, on the intrinsic dynamics of a model pre-BötC neuron over a range of values of tonic synaptic drive conductance (g Tonic ), which tunes excitability. We found that simulated TTX and RZ blockade of I NaP both effectively block intrinsic bursting, either converting bursting dynamics to silence or to tonic spiking, depending in part on the value of g Tonic . Moreover, the effects of these two mechanisms on the shape of the bursting regime and on the frequency and duration of bursts in the appropriate 2D parameter space are nearly indistinguishable (Fig 3).
Since simulated TTX and RZ blockade have remarkably similar effects on I NaP and both effectively abolish intrinsic bursting, it is tempting to conclude that although their mechanisms of action are distinct, TTX and RZ are functionally equivalent I NaP blockers. In the next subsection, we will demonstrate a mathematical similarity in their effects on bursting. In the following subsection, however, we will show how a functional difference in their effects can nonetheless arise.

Fast-slow decomposition analysis
Next, to understand how TTX and RZ affect or abolish intrinsic busting, we used fast-slow decomposition analysis (see [21] for review). This method separates the full system into fast and slow subsystems, with the latter represented by h NaP for the pre-BötC model. Information about the dynamics of the full system can then be inferred from the geometry of the projections of the V m -and h NaP -nullclines into the (V m , h NaP ) phase space. The V m -and h NaP -nullclines are defined from dV m /dt = 0 and dh NaP /dt = 0, respectively. The solution to the latter equation is easy to write down as The V m -nullcline is determined numerically to form a manifold consisting of three branches that connect at folds or "knees", forming an N-shape in the (V m , h NaP )-plane. Intersections between the V m -and h NaP -nullclines represent equilibrium points of the full system. In the bursting regime, the h NaP -nullcline intersects the middle branch of the V m -nullcline and forms an unstable equilibrium point that is surrounded by a stable limit cycle.
Simulated TTX and RZ blockade have distinct effects on the V m -and h NaP -nullclines. Simulated TTX blockade changes the shape of the V m -nullcline by moving the left and right knees to larger and smaller values of h NaP , respectively, although the right knee is not relevant to the behaviors under study. The former effect corresponds to the increased deinactivation of I NaP needed for the neuron to activate after g NaP has been decreased. As g NaP is lowered, the model remains in the bursting regime until the part of the V m -nullcline near the left knee intersects the h NaP -nullcline. The SNIC (saddle-node on an invariant circle) bifurcation that results creates a stable equilibrium point that corresponds to the stabilization of a constant, hyperpolarized membrane potential, marking a transition from bursting, through a decrease in burst frequency, to silence ( Fig 4B, gray dots). In contrast, simulated RZ blockade only affects the h NaP -nullcline, inducing a shift in the leftward, hyperpolarizing direction. In this case, as h NaP is diminished, bursting continues until the h NaP -nullcline intersects the V m -nullcline near the left knee and forms a stable equilibrium point. Interestingly, although RZ blockade affects a different nullcline from TTX, the transition from bursting to quiescence still occurs via a SNIC bifurcation ( Fig 4C, gray dots), which helps explain the similarity in burst properties in

Transient hyperpolarization elicits rebound bursting after I NaP blockade
Although both forms of I NaP blockade convert the model neuron's intrinsic dynamics from bursting to silence via a SNIC bifurcation, we next considered whether transient bursting could be elicited by perturbations. As can be seen in Fig 4A, a burst occurs when h NaP rises to a value above that of the left knee of the V m -nullcline, call it h LK NaP , such that the model's trajectory can transition to high voltages. The position of the entire V m -nullcline, including its left knee, depends not only on the extent of I NaP blockade but also on the level of input to the neuron. Indeed, hyperpolarizing inputs shift the V m -nullcline to more hyperpolarized V m values. As a result, its left branch equilibrium point ends up more hyperpolarized in V m and, due to the voltage-dependence of I NaP , at a more deinactivated, larger h NaP value, call it h EQ NaP . The left knee of the V m -nullcline also rises to a new value,h LK NaP > h LK NaP . In theory, as the model settles toward the equilibrium h EQ NaP , although h NaP remains belowh LK NaP , h NaP could exceed the original left knee height, h LK NaP . If this happens, then release from hyperpolarization will result in a single transient rebound burst. For this result to occur, three conditions need to be met: (1) h LK NaP must be less than 1, so that it is physically realizable, (2) the hyperpolarization must be Pharmacological mechanism and network dynamics sufficiently strong that h EQ NaP > h LK NaP , and (3) the hyperpolarization must be maintained long enough. Fig 5 characterizes the magnitude and duration of hyperpolarization needed to elicit rebound bursting as a function of simulated TTX and RZ blockade of I NaP in our model. We found that with simulated TTX blockade of I NaP rebound bursting is still possible over a narrow range of g NaP , g NaP 2 (1.6 − 3.55 nS). The magnitude and duration of hyperpolarizing input needed to elicit rebound bursting quickly becomes biologically unrealistic, however. In contrast, there are always input magnitudes and durations that yield rebound bursting following any level of simulated RZ blockade, because RZ blockade does not affect h LK NaP and thus condition (1), h LK NaP < 1, holds for all levels of Δh 1/2 . Post-inhibitory rebound bursting has been observed in the intact respiratory network and may play an important role in regulating respiratory rhythms [19,22]. In in vivo conditions, the pre-BötC receives strong inhibition during the interburst interval from other nuclei involved in respiratory rhythm and pattern formation [13,[17][18][19][20]. During this interval the transient hyperpolarization is approximately 10 mV in magnitude and 2 s in duration [13,17]. Our simulations predict that "complete" blockade of I NaP , defined by complete loss of sustained intrinsic bursting (Fig 3), will stop this inhibitory input from inducing a burst in pre-BötC neurons if the block is applied via TTX (Fig 6). This loss occurs because after TTX application, the hyperpolarization does not achieve condition (2), h EQ NaP > h LK NaP . If the block is induced with RZ, however, then the trajectory reaches a position in the (V m , h NaP ) plane where condition (2) does hold, such that subsequent removal of hyperpolarization will result in a burst. The latency (*275 ms) of post-inhibitory rebound bursting in our model is consistent with experimental observed latencies (*300 ms) [22]. Thus, these simulations predict that in the presence of repetitive cycles of inhibition associated with ongoing respiratory rhythms in vivo, pre-BötC neurons will still generate bursts supported by I NaP even after RZ application.
The pre-BötC pre-I network Next, to investigate simulated TTX and RZ application in the pre-BötC network, we constructed a heterogeneous population of 50 model pre-BötC neurons coupled though all-to-all excitatory synapses; see Materials and methods for a full model description. This network is often referred to as the pre-inspiratory/inspiratory (pre-I/I) population and is thought to drive the fictive inspiratory rhythm seen in in vitro slice preparations [13,23,24]. For consistency with experimental observations, parameters were chosen such that approximately 30% of neurons in the synaptically uncoupled network, referred to as "pacemaker" neurons, remain rhythmically active [25]. All other neurons are referred to as "non-pacemaker" neurons. In this network, the type of each neuron (pacemaker vs non-pacemaker) is determined by the neuron's excitability, set by its values of the tonic synaptic conductance g Tonic and the leak reversal potential (E Leak ), the latter of which is normally distributed in order to introduce heterogeneity in the population. I NaP is included in both pacemaker and non-pacemaker neurons and therefore both neuronal types are affected by simulated TTX and RZ blockade. For simplicity, we also tuned the synaptically coupled network such that all neurons, pacemakers and non-pacemakers, are recruited into network oscillations, even though non-pacemakers cannot produce oscillations on their own (Fig 7).

Uniform TTX and RZ block of I NaP in the pre-BötC pre-I network
The synaptic coupling between neurons within the pre-BötC network does not alter the mechanisms of action of TTX and RZ on voltage-gated sodium channels described above. We found that increasing the degree of uniform I NaP block by simulated TTX or RZ results in a progressive reduction in network frequency and a slight reduction in network amplitude, followed by an abrupt cessation of network oscillations (Fig 8A & 8B). With RZ, however, additional off-target effects need to be considered [10,26]. Specifically, at the same concentrations (0 − 20 μM) used to block I NaP in respiratory circuits [9,12,13], RZ also inhibits glutamatergic Pharmacological mechanism and network dynamics excitatory synaptic transmission [11,27]. The magnitude of RZ-mediated attenuation of glutamatergic transmission and its effect on network oscillations within the pre-BötC are not known. Therefore, to understand this off-target effect and allow us to deduce its importance in previous experimental findings (see Fig 9 and related text), I SynE was blocked in a separate set of simulations by systematically reducing the excitatory synaptic conductance parameter (W MaxE ). In contrast to the I NaP effects, we found that simulating the progressive block of synaptic excitation by RZ results in a slight increase in network frequency and a large reduction in network amplitude ( Fig 8C). These outcomes agree with the results from experimental block of excitatory synapses in the pre-BötC [24].
Next, we compared simulated uniform TTX and RZ application with experimental data adapted from [9]. This experimental dataset characterizes the dose-dependent effects of  bilateral microinfusion of TTX and RZ within the pre-BötC of neonatal rat brainstem slices in vitro on the amplitude and frequency of R XII motor output as a function of time. For the comparison of experimental and simulated data, we plotted amplitude versus frequency in order to eliminate the unknown temporal dynamics of TTX and RZ application. Specifically, the experimental data points represent the network frequency and amplitude plotted at successive 1 minute intervals after TTX or RZ application. The corresponding data points plotted from simulated TTX and RZ application represent network frequency and amplitude for increasing levels of blockade. Consistent with the experimental data, we found that simulated TTX application in the pre-BötC progressively reduces the frequency without affecting the amplitude of Experimental data was adapted from [9] and shows the progressive change in amplitude and frequency of network oscillations (monitored by integrated hypoglossal nerve activity) relative to baseline following bilateral microinfusion of TTX or RZ at different concentrations into the pre-BötC. The experimental data points represent the network frequency and amplitude plotted at successive 1 minute intervals after TTX or RZ application. The corresponding data points from simulated TTX and RZ application represent network frequency and amplitude for increasing levels of blockade. Simulated RZ application affects I NaP and excitatory synapses whereas TTX only affects I NaP . In the simulations, the relevant values at the end points (where network oscillations stop) are as follows: (A) g NaP = 3.52; (B) (top trace) Δh 1/2 = −5.0 mV, W maxE = 0.0255 nS, (bottom trace) Δh 1/2 = −4.5 mV, W maxE = 0.0224 nS. Notice that TTX only affects frequency, whereas RZ affects frequency and amplitude. (a1, b1, b2) Effect of simulated TTX and RZ (tunings 1 & 2) blockade on the peak I NaP , peak I Syn , and the I NaP inactivation threshold (h TH NaP ) required to initiate bursting. h TH NaP was defined as the maximal value of the mean population h NaP prior to burst initiation. https://doi.org/10.1371/journal.pcbi.1006938.g009 Pharmacological mechanism and network dynamics population oscillations before oscillations abruptly stop ( Fig 9A). With RZ, however, assuming that both I NaP and excitatory synapses are impacted, extrapolation from our separate simulations of these effects predicts that experimental RZ application will progressively reduce both frequency and amplitude of the pre-BötC population oscillations until they eventually stop. Fig  9B demonstrates that these effects do arise both in simulated RZ application, where both I NaP and synapses are affected, and in experimental data [9]. We found that our simulation results matched the reduction in network amplitude and frequency seen with experimental RZ application when W MaxE decayed exponentially with increasing hyperpolarization of h NaP 1=2 (see Materials and methods). To match experimental results of RZ application at 5 μM, W MaxE was maximally reduced by 7% (Fig 9B, Tuning 1), and for RZ at both 10 μM and 20 μM, W MaxE was maximally reduced by 15% (Fig 9B, Tuning 2).
To investigate the mechanisms underlying changes in network amplitude and frequency, we quantified changes in the I NaP threshold (I TH NaP ) required for spike generation and the peak I SynE magnitude during network bursts, and we estimated the h NaP threshold (h TH NaP ) required to initiate bursting as a function of the level of TTX and RZ application in our simulations, tuned to match experimental data from [9] (Fig 9a1, 9b1 and 9b2).
In the model, amplitude is defined as the number of spikes per neuron per 50 ms bin. Therefore, by this definition, changes in network amplitude result from changes in the firing rate of bursting neurons and/or changes in the number of neurons recruited into network bursts. The number of recruited neurons is affected by I SynE and the firing rate of individual neurons during bursting is a function of the total depolarizing current, determined by both I NaP and I SynE . In contrast, network oscillation frequency is determined by the time required for I NaP to recover from inactivation to a sufficient threshold (I TH NaP ) for spike generation following each network burst. To be more precise, network frequency is a function of the time required for h NaP to recover from inactivation up to a level, h TH NaP , such that I TH NaP is reached following each network burst. For a fixed level of synaptic input, threshold I TH NaP is a constant given by I TH NaP ¼ À g NaP � m NaP � h TH NaP � ðV m À E Na Þ (see Eq 5 in Materials and methods). Thus, as long as synaptic input remains constant, h TH NaP increases as g NaP decreases. In the period of h NaP recovery leading up to burst onset, we indeed have the constant synaptic input I SynE = 0, independent of W MaxE , since the pre-BötC neurons activate synchronously in the regime we are considering.
In the case of simulated TTX, while I TH NaP and I SynE did not change, the above reasoning yields an increase in h TH NaP with the progressive decrease in g NaP (Fig 9a1). Note that the increase in the h TH NaP is consistent with Fig 4. Therefore, the reduced network frequency seen with increasing TTX blockade is a direct consequence of an increased h NaP recovery time driven by the increased h TH NaP and the proximity of the nullclines (Fig 4). Because I TH NaP does not vary with g NaP , the level of I NaP expressed in bursting neurons is the same as before the blockade. This invariance of I NaP together with the fact that TTX does not affect I syn imply that the firing rate of bursting neurons and the number of recruited neurons are unchanged. This reasoning explains why network amplitude is unaffected by simulated TTX blockade of I NaP .
In the case of simulated RZ, since g NaP does not vary, h TH NaP is unaffected by progressive blockade. In this scenario, changes in frequency arise instead due to a slower rate of recovery from inactivation caused by the hyperpolarizing shift in the inactivation curve (i.e., the h NaPnullcline), Figs 2 & 4. This shift is accompanied by a reduction of I SynE that follows directly from the reduction of g SynE with simulated RZ. Therefore, changes in network amplitude with simulated RZ arise due to the progressive reduction of I SynE , which results in a decreased firing rate during bursting and de-recruitment of a subset of the neurons.

Non-uniform TTX and RZ block of I NaP in the pre-BötC pre-I network
Blockade of I NaP by experimental application of TTX or RZ is likely to be non-uniform, especially in the case of in vitro bath application where drug penetration depends on passive diffusion. With TTX and RZ bath application, I NaP will be blocked in neurons closest to the surface before neurons at the center of the slice. As mentioned previously, the pre-BötC pre-I population contains pacemaker and non-pacemaker neurons that play different roles in rhythm and pattern formation. The spatial orientation of these two types of neurons within the pre-BötC is unknown. Therefore, to understand the effects of non-uniform drug penetration in this model, we simulated the sequential blockade of I NaP by TTX and RZ in the pre-I population, under each of three different assumptions on the order in which pacemaker and non-pacemaker neurons are affected: (1) pacemaker then non-pacemaker, (2) non-pacemaker then pacemaker, and (3) random order. We found that the effects of TTX and RZ are indistinguishable for all three cases. If pacemaker neurons are affected first, then I NaP block via either mechanism results in a progressive reduction in network frequency with no effect on amplitude until oscillations eventually stop ( Fig 10A). That is, with loss of I NaP in pacemakers, burst initiation is slowed, but once a burst occurs, all neurons are recruited. In contrast, if non-pacemaker neurons are affected first, then I NaP block results in a progressive reduction in network amplitude and an increase in frequency before oscillations eventually stop ( Fig 10B). The loss of I NaP in non-pacemakers can prevent their recruitment, while the involvement of fewer neurons in each burst results in shorter bursts that can recur more frequently. Finally, if I NaP is blocked in pacemaker and non-pacemaker neurons in a random order, then I NaP block results in a progressive reduction in network amplitude and frequency before oscillations eventually stop. In all cases, network oscillations stop as soon as I NaP is completely blocked in all pacemaker neurons.

Simulated TTX and RZ block of I NAP in the intact respiratory network
Finally, to investigate the effects of simulated TTX and RZ block of I NaP in the intact respiratory network, we reconstructed a multi-population network model that is thought to represent the core mammalian respiratory central pattern generator located in the BötC and pre-BötC [13,17,[28][29][30][31]. The simulated intact network is composed of four subpopulations each consisting of 50 neurons. The subpopulations are characterized by the timing of their activity relative to the phases of inspiration and expiration as follows: post inspiratory (post-I), augmenting expiratory (aug-E), pre-I/I as simulated in the pre-BötC network considered thus far in this work, and early-inspiratory (early-I). This network produces a three-phase rhythm with inspiratory (I), post-inspiratory (pI), and stage-2 expiratory (E2) phases, which is similar to respiratory rhythms seen in vivo (Fig 11). For a full model description, see Materials and methods.
In the intact respiratory network, the respective contributions of I NaP and inhibitory network interactions to rhythm generation vary depending on the excitability state of the pre-I population [29]. Therefore, before simulating TTX and RZ blockade of I NaP in the intact network, we first characterized the dependence of pre-I population bursting on I NaP and synaptic inhibition in the intact network as a function of g Tonic (Fig 12A). This was accomplished by comparing the dynamics of the pre-I population (silent, bursting, tonic), embedded in the network, between baseline conditions and the extreme cases where I NaP = 0 or I SynE = 0. We found that under baseline conditions, bursting in the pre-I population is extremely robust and oscillations continue over all tested values of g Tonic , namely 0.3 − 1.0 nS (cf. [32]). Complete block of inhibitory synaptic currents reveals that the pre-I population bursting is dependent on inhibitory network interactions for g Tonic 2 (0.40, 1), indicated by a transition of the pre-I population from bursting, under baseline conditions, to a tonic spiking mode of activity when I SynE = 0. In contrast, complete blockade of I NaP reveals a regime of I NaP -dependent bursting for g Tonic = 0.3 − 0.486, indicated by a transition of the pre-I population from a bursting to a silent mode of activity when I NaP = 0. Importantly, these two regimes overlap (for g Tonic 2 (0.40, 0.486)), which indicates a region where both I NaP and inhibitory synaptic interactions are critical for rhythm generation/bursting in the pre-I population (Fig 12A).
Next we characterized network parameters affecting the magnitude/contribution of I NaP during pre-I population bursting, which, as in the pre-I network, modulate I NaP inactivation dynamics. As with the pre-I network, there is a threshold level of I NaP needed for pre-BötC burst onset, which is realized when h NaP deinactivates to a corresponding threshold level, h Th NaP . One new feature in this case, however, is that h Th NaP depends both on g Tonic and on the level of  Pharmacological mechanism and network dynamics synaptic inhibition to the pre-BötC population, since this parameter and input affect the position of the V m -nullcline for the pre-BötC neuron (analogously to the impact of applied inhibitory current on a single pre-BötC neuron shown in Fig 6B). This dependence holds in both the I NaP -dependent and the network-dependent burst regimes (Fig 12). Therefore, in the pre-I population, we characterized: (1) the maximal hyperpolarization during the fictive expiratory phase, which impacts the rate of deinactivation of h NaP , (2) the dynamic range of h NaP over a complete respiratory cycle, and (3) the peak I NaP during bursting, as a function of the strength of inhibition from the post-I population for three levels of pre-I population tonic excitation set by g Tonic (Fig 12B, 12C and 12D).
Increasing g Tonic naturally decreases the hyperpolarization of the pre-I neurons in the expiratory phase (Fig 12B). Due to its effects on the V m -nullclines for the pre-I neurons, it also lowers I Th NaP and hence h Th NaP . These changes lower the operating range of h NaP (Fig 12C) and correspondingly result in a lower level of I NaP while the pre-I population is active (Fig 12D). In contrast, increasing the strength of inhibition from the post-I population to the pre-I neurons hyperpolarizes them more ( Fig 12B) and increases I Th NaP and hence h Th NaP through the effect of inhibition on the pre-I V m -nullcline. These changes raise the operating range of h NaP to higher levels ( Fig 12C) and correspondingly allow I NaP to reach a higher level. Note that the contribution of I NaP is small when g Tonic is relatively strong and post-I inhibition is weak, whereas the contribution of I NaP is large when the relative strengths of g Tonic and post-I inhibition are reversed. Interestingly, under the latter conditions, the magnitude of I NaP can be larger in the intact network than in the isolated pre-I population due to the post-I population inhibition, which can enhance I NaP deinactivation.
Finally, we simulated TTX and RZ blockade of I NaP in the intact respiratory network ( Fig  13). The effects of simulated I NaP blockade on the pre-I population dynamics are expected to depend on the value of g Tonic and the strength of post-I inhibition within this population, due to the effects illustrated in Fig 12. Therefore, we characterized the relative change in the pre-I population amplitude and frequency as a function of g Tonic and post-I inhibition for 'complete' blockade of I NaP by simulated TTX and RZ. Blockade of I NaP was considered complete for simulated TTX when g NaP = 0 nS, and for simulated RZ when Δh 1/2 = −15 mV (see Fig 2).
In the case of simulated TTX, complete blockade of I NaP abolishes pre-I population oscillations for g Tonic < 0.486 nS across all levels of post-I inhibition considered (Fig 13A), as expected from the boundary between I NaP -and network-dependent bursting regions in Fig  12A. In contrast, for g Tonic > 0.486 nS pre-I population oscillations continue after complete blockade of I NaP . In this regime, network amplitude is generally reduced without affecting frequency, except for a frequency reduction for g Tonic � 0.486 nS. The extent of the reduction in amplitude depends on g Tonic and the strength of post-I inhibition. The reduction in amplitude is largest when g Tonic is weak and post-I inhibition is strong.
In the case of simulated RZ, the range of g Tonic values where complete blockade of I NaP abolishes pre-I population bursting is greatly reduced (Fig 13B) compared to simulated TTX. In general, stronger post-I inhibition decreases the range of g Tonic values where I NaP block abolishes pre-I population bursting. This relationship is due to the increased recovery of I NaP from inactivation caused by the strengthening of post-I inhibition and the enhanced expiratory phase hyperpolarization of pre-I that it induces. This mechanism is illustrated in Fig 6. In the region of parameter space where pre-I population bursting is not abolished, network amplitude is generally reduced by RZ blockade without affecting frequency. An exception occurs close to the boundary where bursting is almost abolished; in this case RZ blockade reduces both amplitude and frequency. Moreover, the relative decrease in amplitude is largest when g Tonic � 0.4 nS and post-I inhibition is strong (≳ 0.5 nS), Fig 13. Importantly, in these simulations, complete blockade of I NaP by RZ fails to stop the respiratory rhythm under network conditions where pre-I population bursting is dependent on I NaP (Fig 12A).
Finally, we compared the effects of simulated RZ blockade of I NaP in the intact respiratory network with experimental data from [13] (Fig 13C). This experimental data set characterizes the steady-state dose-dependent effects of RZ on burst amplitude and frequency of the intact respiratory network output measured from integrated phrenic nerve activity. In these Pharmacological mechanism and network dynamics experiments, the maximal 20 μM RZ concentration, which is assumed to completely block I NaP , resulted in a * 35% reduction in network amplitude and no significant change in frequency. In order to match these data in our simulations, we found that post-I inhibition must be relatively strong (>5 nS) and g Tonic between 0.375 nS and 0.425 nS (Fig 13B and 13C). Interestingly, this range of g Tonic values overlaps with the regimes where pre-I population bursting is exclusively I NaP -dependent and where bursting is dependent on both I NaP and inhibitory network interactions (Fig 12A). Consequently, under these network conditions, our model predicts that selective and complete block of I NaP by experimental application of TTX at low concentrations (� 20 nM) will abolish the respiratory rhythm in the intact network (Fig 13D).

Discussion
Understanding how pharmacological mechanisms and network dynamics affect I NaP blockade by TTX and RZ is critical for interpreting experimental data and its implications for understanding the underlying mechanisms of respiratory rhythm and pattern formation. Therefore, in this computational study, we characterized the effects of TTX and RZ blockade of I NaP on respiratory network dynamics by simulating their distinct pharmacological mechanisms of action in established models of respiratory neurons and neurocircuitry that represent in vitro and in vivo mouse/rat preparations. To summarize, we show that in simulated pre-BötC respiratory neurons under conditions representing in vitro slice preparations, TTX and RZ both effectively block I NaP and abolish intrinsic bursting (Figs 2 & 3). Given these findings, it is tempting to conclude that TTX and RZ application will induce similar effects on respiratory dynamics under all experimental conditions. Interestingly, however, we found that after simulated blockade of RZ, but not TTX, I NaP can be reactivated by transient hyperpolarization, due to differences between these drugs' specific pharmacological mechanisms of action (Fig 6). This effect becomes critical in the intact respiratory network where these neurons receive strong, transient hyperpolarizing inhibition from post-I neurons during the expiratory phase of respiration. Correspondingly, our simulations of I NaP blockade in the intact respiratory network predict that experimental application of RZ, but not TTX, will fail to effectively block I NaP and respiratory rhythmicity.

Insights into the role of I NaP in the intact respiratory network
This computational study provides novel insight into the role of I NaP in respiratory rhythmogenesis and pattern formation. In the intact respiratory network, the role of I NaP is not well understood. Current thinking is that in the intact network I NaP is largely inactivated and hence not essential for ongoing rhythm generation and pattern formation [13,17]. This idea is based on computational modeling and the experimental observation that RZ application in the intact network fails to stop the respiratory rhythm, suggesting that active recruitment of I NaP may not be essential to rhythmogenesis under these conditions [13].
This interpretation does not consider the dynamic interaction between network inhibition, tonic excitation and the voltage-dependent inactivation dynamics of I NaP , however. Moreover, this interpretation assumes that RZ effectively blocks I NaP and overlooks the experimental observation that RZ application significantly reduces the amplitude of the respiratory rhythm generated by the intact network [13]. Indeed, the latter observation is not consistent with the idea that I NaP is largely inactivated, as blocking an inactivated current should have no effect on rhythm characteristics. Our simulations support an alternative view. Specifically, to generate comparable reductions in network amplitude from simulated RZ blockade based on a shift in the I NaP inactivation curve and from simulated TTX blockade based on a reduction in I NaP conductance, I NaP must strongly contribute to pre-I population bursting. Importantly, this contribution only occurs if, during periods when pre-I neurons are not spiking, I NaP strongly deinactivates and h NaP levels are higher in the intact network than in the isolated pre-I population (Figs 12 & 13).
An important finding of this study is the conclusion that I NaP plays a critical role in respiratory rhythmogenesis, which is indicated by the prediction that complete blockade of I NaP will stop respiratory rhythm generation in the intact network. Importantly, this finding comes to light only by considering the distinct pharmacological mechanisms of TTX and RZ in the context of the interaction between pre-I population excitability, inhibitory network interactions, and the associated dynamics of I NaP inactivation.
We showed that the dependence of pre-I bursting on I NaP in the intact network is a function of the excitability of the pre-I population (Fig 12A). When inputs or neuromodulation sufficiently lowers the excitability of pre-I neurons, I NaP becomes a necessity for burst dynamics, whereas with heightened excitability, bursting is dependent on inhibitory network interactions. While these mechanisms can act separately, there is a range of excitabilities for which pre-I bursting depends on both I NaP and inhibitory network interactions.
In the intact network, the pre-I population receives transient inhibition from post-I neurons during the expiratory phase of respiration [13,[17][18][19][20], which, our results show, may compromise the ability of RZ to effectively block I NaP . Our RZ simulations, tuned to match experimental data (specifically, a large reduction in pre-I output amplitude and no effect on frequency resulting from RZ application), suggest that the excitability of the pre-I population in the intact network is in a regime where rhythm generation depends on both I NaP and inhibitory network interactions (Fig 12A & 13). In these simulations, since RZ blockade of I NaP fails to stop the rhythm and pre-I bursting is I NaP -dependent, this indicates that RZ fails to completely block I NaP . In contrast, our simulations predict that under these network conditions a complete blockade of I NaP by TTX will stop the respiratory rhythm. Importantly, at low concentrations (� 20 nM), TTX has been shown to selectively block I NaP without affecting the action potential generating fast Na + current [9]. Therefore, this prediction can be experimentally tested via bilateral microinfusion of TTX at low concentration into the pre-BötC in the intact preparation. If confirmed, this finding would illustrate the critical role of I NaP in rhythm generation within intact respiratory circuits.
One caveat here is that the extent to which rhythm generation in the intact network relies on I NaP and on inhibitory network interactions is highly dependent on the excitability of the pre-I population (Fig 12). Our simulations, tuned to match the results of [13], suggest that the excitability in the pre-I population must be close to the border between a regime where pre-I bursting depends exclusively on I NaP and a regime where pre-I bursting depends on both I NaP and inhibitory network interactions (Figs 12 and 13). Therefore, with reasonable levels of variability, complete block of synaptic inhibition under these network conditions will inconsistently stop pre-I population oscillations and in instances where oscillations are stopped, the pre-I population will transition to a tonic mode of activity. Indeed, the prediction of proximity to the border is consistent with recent experimental data [18], which showed that simultaneous block/attenuation of GABAergic and glycinergic synaptic inhibition within the pre-BötC failed to abolish rhythmic phrenic nerve output in 83% experiments, and in experiments where oscillations were stopped, phrenic nerve activity culminated in tonic activity. Conditions that increase the excitability state of the pre-I population may transition the intact network into a regime where rhythm generation is no longer I NaP -dependent, and complete I NaP blockade under these conditions is not predicted to stop rhythmic oscillations. Hypoxia and hypercapnia are perturbations that are likely to affect pre-I excitability and consequently the role of I NaP in the intact network.
It is important to note that the model used in the current study is the same as that considered in a past work that came to some different conclusions [13]. Specifically, the authors of that work suggested that I NaP is largely inactivated and hence is not critical for rhythm generation in the intact network. This previous study, however, did not systematically consider the interaction between pre-I population excitability, inhibitory network interactions, and the resulting dynamics of I NaP inhibition. It also did not consider the distinct pharmacological mechanism of action of RZ and I NaP block was simulated by reducing g NaP . With the parameter set used for simulations representing the intact network in [13], the pre-I population excitability was set at g tonic = 0.98 nS and the strength of post-I to pre-I inhibition was 0.225 nS. In [13] and our simulations, I NaP is largely inactivated under these conditions. Consequently, under these parameters, simulated TTX or RZ block of I NaP cannot capture the large reduction in amplitude seen with experimental RZ application in the intact network. Therefore, conclusions drawn from these simulations about the excitability state of the pre-I population, I NaP inactivation dynamics, and the role of I NaP in rhythm generation in the intact respiratory network may not accurately represent the underlying mechanisms, dynamics, and conditions in these experimental preparations.

RZ-dependent reduction in excitatory synaptic transmission
Characterizing the off-target effects of RZ in respiratory circuits may be critical for understanding how RZ application impacts pre-I network dynamics [26]. Comparison of experimental and simulated data in this study supports the idea that RZ (but not TTX) application not only alters I NaP but also attenuates excitatory synaptic transmission. In regions outside of the pre-BötC, RZ has been shown to block excitatory synaptic transmission at doses comparable to those used in respiratory circuits (0 − 20 μM) [11,27]. The effects of RZ on synaptic transmission or other off-target effects have not been characterized in the pre-BötC, however. Experimentally, bilateral microinfusion of TTX or RZ into the pre-BötC stops the fictive respiratory rhythm in in vitro slice preparations. Analysis of the time course of these drugs reveals that before rhythm termination, RZ reduces the frequency and amplitude of integrated hypoglossal nerve motor output, whereas TTX only reduces frequency [9].
In our simulations of the isolated pre-I network, selective blockade of I NaP by either the pharmacological mechanism of TTX or RZ reduces network frequency without affecting amplitude (Fig 8A and 8B). In contrast, blockade of excitatory synaptic transmission selectively reduces network amplitude with little effect on frequency (Fig 8C), which is consistent with experimental data [24]. Therefore, to account for the reduction in amplitude seen with experimental RZ application, our results suggest that RZ must reduce the strength of the excitatory synaptic transmission in addition to modulating I NaP (Fig 9). The prediction that RZ attenuates I SynE is testable by isolating the synaptic current in vitro using voltage-clamp recordings in conjunction with RZ and TTX application.
Importantly, although RZ blocks both I NaP and I SynE in our simulations of the isolated pre-I network, rhythm generation in this system is abolished due to the reduction of I NaP , not I SynE . Partial reduction of I SynE alone is not sufficient to stop rhythm generation in our simulations (Fig 8C). These results support the hypothesis that I NaP is critical for rhythm generation in the pre-BötC in in vitro slice preparations.
It is important to mention that additional off-target effects of RZ have been reported in neurons outside of respiratory circuits such as: potentiation of calcium-dependent K + current, inhibition of fast Na + current, inhibition of voltage-gated Ca 2+ current, inhibition of voltagegated K + current, and inhibition of the glutamate receptor N-methyl-D-aspartate receptor (NMDA) [10,26]. Detailed consideration of these additional off-target RZ effects are beyond the scope of this study. Many of these off-target effects occur at relatively high RZ concentrations and are therefore not likely to affect the results of this study.
Note that our simulations of I NaP block in the isolated pre-I population are specifically tuned to match data where TTX or RZ application is delivered via bilateral microinfusion into the I NaP , as opposed to bath application [9]. This choice was made because bilateral microinfusion is likely a more targeted/localized approach and avoids the issues of non-uniform and offtarget drug effects that may arise with bath application.

Non-uniform pharmacological blockade in the pre-BötC pre-I population
The role of I NaP in respiratory rhythm generation in the isolated pre-I population is also unclear and a highly debated topic within the field. This debate is fueled by the diverse and seemingly contradictory effects of I NaP block by bath application of TTX or RZ in in vitro slice preparations. For example, [9] found that in vitro bath application of either TTX or RZ reliably decreased the frequency and to a lesser extent the amplitude of network bursting before oscillations abruptly stopped. In contrast, [25,33] found that bath application of RZ at comparable concentrations resulted in a large reduction in the network amplitude but had no effect on or even increased frequency and failed to stop network oscillations. Consideration of non-uniform I NaP block may provide an explanation for these seemingly contradictory findings.
With bath application of TTX or RZ, block of I NaP is unlikely to be uniform since drug penetration is dependent on passive diffusion. Therefore, with bath application, I NaP block will affect neurons close to the surface of the slice first and progress towards neurons in the center. In Fig 10 we show that the result of progressive I NaP block across the isolated pre-I population is highly dependent on the dynamics of the neurons (pacemaker vs non-pacemaker) affected. Our simulations can explain the diversity of experimental results if we consider the effects of slice thickness and assume that pacemaker neurons are preferentially located near the center of the larger population of pre-I neurons; see Pharmacological mechanism and network dynamics experimental preparations. For example 630 − 690 μm slices were used in [25] and 250 − 350 μm slices are used in [9]. The assumption that pacemaker neurons are preferentially located near the center of the larger population of inspiratory neurons is supported by a recent study [34]. This study found that the size of the active inspiratory (pre-I) network dynamically expands/contracts around a core rhythmogenetic center by recruiting otherwise inactive (nonpacemaker) neurons into pre-I population oscillations along a rostrocaudal axis in response to behavioral or metabolic challenges to breathing.
In a thin slice, both pacemaker and non-pacemaker neurons are located close to the surface and an applied drug should fully diffuse through the slice. In this case, as TTX or RZ diffuses in, I NaP will be blocked in both pacemaker and non-pacemaker neurons, which is predicted to affect both amplitude and frequency until oscillations eventually stop (Fig 10C). In contrast, with a thick slice, TTX and RZ may not diffuse all the way to the center and the neurons affected first would primarily be the peripheral non-pacemaker neurons. With I NaP blocked predominantly in non-pacemaker neurons, our results predict a decrease in the amplitude of oscillations (Fig 10B). If the drug diffuses deeper into the slice, both pacemaker and non-pacemaker neurons will be impacted, which will further decrease amplitude and start to decrease frequency ( Fig 10C). Therefore, our model predicts that in the case of a thick slice, the lack of frequency effects and the failure to stop network oscillations derive from incomplete/inadequate drug penetration into the slice. These predictions are consistent the experimental results of [25] and [9], which use thick (630 − 690 μm) and thin slice preparations (250 − 350 μm), respectively. It is unknown how far TTX or RZ diffuse into neuronal tissue. However, a straightforward prediction of incomplete drug penetration is the observation of pacemaker neurons that appear insensitive to I NaP block. Consistent with this prediction, [25] found pacemaker neurons that are insensitive both to RZ application and to Cd 2+ application that blocks Ca 2+ dynamics thought to underlie rhythmic properties in some pacemaker neurons. An experimental preparation that contains the same population of neurons and likely avoids the issue of non-uniform/incomplete I NaP block is an in situ rat brain stem-spinal cord preparation where regions rostal to the pre-BötC are removed and RZ is delivered by arterial perfusion. In this preparation, application of RZ at concentrations greater than 10 μM consistently stopped rhythmic oscillations of the pre-I population [13]. In addition to providing a mechanistic explanation to seemingly contradictory experimental data, these simulations demonstrate the importance of considering the spatial/temporal dynamics of pharmacological diffusion into thick and thin in vitro slice preparations when interpreting data from pharmacological blocking studies.

Limitations of this study
In the pre-BötC, coupling of a calcium-activated non-selective cation current (I CAN ) and intracellular calcium transients is proposed to represent a biophysical mechanism underlying rhythmogenesis at the cellular-and/or network-level [25,[35][36][37][38][39][40][41]. For simplicity, however, explicit representation of I CAN and intracellular calcium dynamics were omitted in this study. This choice can be justified by considering the effects of pharmacological block of I CAN in conjunction with recent data-driven computational work. In in vitro slice preparations containing the pre-BötC, pharmacological blockade of I CAN results in a large reduction in the amplitude of network oscillation with no or minor perturbation of frequency, suggesting that the primary role of I CAN in these circuits is amplitude generation/regulation [25,42]. A recent data-driven computational study that included I NaP , I CAN , and intracellular calcium dynamics [43] found that in order to match experimental data, I CAN activation must be strongly coupled to synaptically triggered calcium transients, as suggested by [37,[39][40][41]44]. Consequently, I CAN acts as a mechanism that amplifies excitatory synaptic inputs within the pre-BötC. As such, I CAN can be treated as an excitatory post synaptic current that, in our model, could correspond to a portion of the total excitatory current I SynE . Consistent with this idea, blockade of I SynE in our model of the isolated pre-BötC is consistent with the effects seen with experimental blockade of I CAN (Fig 8C).
Another issue that we do not explore in this work is the contribution of additional burst termination mechanisms. At the neuronal and population level, mechanisms of burst termination are critical for the generation of ongoing rhythmic oscillations. In the current model, burst termination is exclusively dependent on the inactivation of I NaP . However, in respiratory neurons and circuits, additional mechanisms of burst termination have been proposed, such as slowly activating potassium channels [14], inositol triphosphate (IP3) effects [38,45,46], Na + /K + ATPase electrogenic pumps [46,47], and synaptic depression [40,48], for example. Inclusion of additional burst terminating mechanisms in I NaP -dependent rhythmogenic models of pre-BötC neurons increases the dynamic range where bursting occurs, making rhythm generation more robust [46], and may change the effects of I NaP block on neuronal activity in these models. For example, inclusion of a Na + /K + ATPase electrogenic pump in [46], in addition to I NaP , dramatically increases the dynamic range where I NaP block transitions the neuron from a bursting to a tonic mode of activity, as opposed to a transition to a silent model of activity. Although theoretically interesting, additional burst termination mechanisms were not necessary to match a wide range of experimental data in our simulations and therefore were not included in this study.
Several additional properties may modulate pre-BötC oscillations but are not considered in this work. Local inhibition within the pre-BötC, for example, may play a role in pre-I burst synchrony, variability and patterning [20], but is not critical for intrinsic pre-BötC oscillations [49]. Similarly, spatial variability in network organization and synaptic connection probability as well as stochasticity and neuromodulation may play important roles in shaping pattern formation in respiratory circuits [20,39,[50][51][52][53], but consideration of these additional features is beyond the scope of this study and is left for future work.

Broader implications
I NaP is not unique to respiratory circuits. Experimental and computational results have suggested an active role for I NaP in shaping activity patterns of putative locomotor CPG neurons [54][55][56] and in contributing to locomotor rhythm generation [55,57,58]. Similar roles are supported for I NaP in generating the neural rhythm underlying mastication [59,60]. I NaP is also expressed and contributes to neural activity elsewhere in the central nervous system [61,62]. RZ and TTX have been used in experimental investigations of these neurons and their activity patterns. Additionally, RZ has been considered as a neuroprotective or anticonvulsant agent [8] and is actively investigated for its therapeutic potential in the treatment of multiple diseases such as obsessive-compulsive disorder (OCD) [63][64][65], anxiety [66,67], depression [68,69], multiple sclerosis (MS) [70], and amyotrophic lateral sclerosis (ALS) [71]. Therefore, understanding RZ's pharmacological mechanisms of action may be important for understanding its beneficial effects in these disease states.
Additionally, the specific pharmacological mechanisms of action considered in this study are not unique to TTX and RZ. Obstruction of the ion channel pore and shifts in voltagedependent activation/inactivation dynamics are common pharmacological mechanisms of action [1,2]. Similarly, spatial/temporal dynamics of drug diffusion in thick and thin in vitro slice preparations are not unique to the medulla. Therefore, the findings of this study are broadly relevant when interpreting and simulating data from pharmacological blocking studies.

Conclusion
This computational study, in conjunction with experimental data, has illustrated the general importance of considering the pharmacological mechanism(s) of action when interpreting and simulating data from pharmacological blockade studies. In the case of respiratory circuits, we show that (1) RZ may fail to effectively block I NaP in the intact network due the specific pharmacological mechanism of action and transient inhibitory network interactions, (2) pre-I bursting in the intact network is likely dependent on both inhibitory network interactions and I NaP , and (3) experimental TTX application in the intact network is predicted to terminate respiratory rhythmicity. These findings suggest that I NaP is less inactivated than has been previously proposed and that it plays a critical role in respiratory rhythm generation under in vivo conditions. If experimentally confirmed, these predictions will advance our understanding of the mechanisms of rhythm generation in brainstem respiratory circuits. These findings also have implications in understanding the role of I NaP in CPGs other than respiratory circuits, in elucidating the effects of RZ across the CNS in the treatment of various conditions, and in interpreting and simulating data from pharmacological blockade studies in general.

Model description
Neurons were simulated with single compartment models incorporating Hodgkin-Huxley style conductances based on previously described models [13,14,46]. The membrane potential V m for each neuron is given by the following differential equation: where C m = 36.0 pF is the cell capacitance, I Na , I K , I Leak , I NaP , I Ca , I KCa , I SynE , and I SynI are the sodium, potassium, leak, persistent sodium, high-voltage activated calcium, calcium-activated potassium, excitatory synaptic and inhibitory synaptic ionic currents, respectively. The currents are defined as follows: I NaP ¼ g NaP � m NaP � h NaP � ðV m À E Na Þ ð5Þ where g i is the maximum conductance, E i is the reversal potential, and m i and h i are gating var- where a 1 ðVÞ ¼ A a � ðV þ B a Þ=ð1 À expðÀ ðV þ B a Þ=k a ÞÞ; ð15Þ The values for the constants A α , A β , B α , B β , k α , and k β are also given in Table 1.
For I KCa , the steady-state activation (m KCa 1 ðCa in Þ) and time constant (t KCa m ðCa in Þ) depend on the intracellular Ca 2+ concentration (Ca in ) as follows: The values of α KCa , β KCa and t m max can be found in Table 1. Ca in is determined by the balance of Ca 2+ influx carried by I Ca and efflux via the Ca 2+ pump. The dynamics of Ca in is described as follows: where α Ca = 2.5 � 10 −5 mM/fC is a conversion factor relating current to rate of change of Ca in , τ pump = 500ms is the time constant for the Ca 2+ pump and Ca min = 5.0 � 10 −6 mM is a minimal baseline calcium concentration. The total synaptic conductances (g i SynE and g i SynI ) of the i th target neuron are described by the following equation: where g i Tonic is the tonic excitatory conductance, W E j;i and W I j;i are the weights of the excitatory and inhibitory synaptic connection from the source neuron j to the target neuron i, H(.) is the Heaviside step function, and t denotes time. τ SynE and τ SynI are exponential decay constants for excitatory and inhibitory synapses. t j,n is the time at which the n th action potential is generated in neuron j and reaches neuron i.

Network construction
The simulations representing the intact respiratory network and the isolated pre-I population were reconstructed from the model description used in [13]. In the intact network each population (pre-I, early-I, aug-E, post-I) consists of 50 model neurons governed by the equations presented in the previous subsection, with intrinsic neuronal parameter values given in Table 2. Note that heterogeneity was introduced by uniformly distributing the parameters E leak within each population as well as the weights of excitatory (W E j;i ) and inhibitory (W I j;i ) synaptic connections; see Tables 2 & 3. The tonic excitatory drive (g Tonic ) within the isolated pre-I population model (meant to represent in vitro slice preparations) was tuned such that the percentage of bursting neurons in the synaptically decoupled network is between 20-30%, consistent with experimental observations [25]. The maximal weight (W MaxE ) of excitatory synapses between neurons was then tuned to achieve complete network synchronization, where all neurons are active during network oscillations. In full network simulations (meant to represent in vivo or in situ brain stem-spinal cord preparations), g Tonic was treated as a variable to match experimental observations [13].

Data analysis and definitions
Data generated from simulations was post-processed in Matlab (Mathworks, Inc.). An action potential was defined to have occurred in a neuron when its membrane potential V m increased through −35mV. Histograms of population activity were calculated as the number of action potentials per 50ms bin per neuron with units of APs/(s � neuron). Population amplitude and frequency were calculated by identifying the peaks and the inverse of the interpeak interval from the population histograms.

Integration methods
All simulations were performed locally on an 8-core Linux-based operating system. Simulation software was custom written in C++. Numerical integration was performed using the firstorder Euler method with a fixed step-size (Δt) of 0.025ms. Pharmacological mechanism and network dynamics