article: 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.
  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 (INaP) 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 INaP inactivation in the hyperpolarizing direction [4–6]. At low concentrations, both TTX (≤ 20 nM) and RZ (≤ 20 μM) have been shown to selectively block INaP over the fast action-potential generating sodium current (INa) [7–10]. A caveat to RZ blockade, however, is that RZ has been shown to inhibit excitatory synaptic transmission at concentrations that affect INaP [10,11]. In respiratory circuits, in vitro blocking studies using TTX or RZ have suggested that INaP 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 INaP may be a necessary component of rhythmogenesis in respiratory circuits [14–16]. This hypothesis, however, has fallen out of favor due to the observation that INaP block by RZ fails to stop respiratory rhythms in intact preparations [13]. This conclusion is dependent on the assumption that RZ effectively blocks INaP in an in vivo setting. This study predicts that after blockade via RZ, INaP can be reactivated by transient hyperpolarizing perturbations due to the specific mechanism of action of RZ. Simulated TTX blockade of INaP does not yield this same effect. In intact preparations, the pre-BötC receives strong inhibition during the interburst interval [13,17–20], which may allow INaP 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 INaP, while complete block of INaP 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 INaP 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.
  Bursting in this model is dependent on INaP, therefore, we first characterized the effects of simulated TTX and RZ blockade on INaP under steady-state conditions (Fig 2). Since TTX directly obstructs the pores of sodium channels, TTX blockade was simulated by reducing the conductance gNaP of INaP. Since TTX only affects conductance, the steady-state activation and inactivation parameters of INaP were not varied (Fig 2A1). We observed that the steady-state current is reduced proportionally to the decrease of gNaP (Fig 2A2). In contrast, RZ shifts voltage-dependent inactivation in the hyperpolarizing direction. Therefore, RZ blockade was simulated by decreasing the half-inactivation parameter h1/2 (Fig 2B1). We found that, similar to TTX, simulated RZ application effectively blocks INaP in steady-state conditions (Fig 2B2). Interestingly, although TTX and RZ blockade of INaP work through distinct mechanisms, the effects on the peak steady-state current are remarkably similar (Fig 2C). The activity pattern exhibited by a neuron capable of INaP-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 INaP, represented by reduction of gNaP or h1/2, respectively, on the intrinsic dynamics of a model pre-BötC neuron over a range of values of tonic synaptic drive conductance (gTonic), which tunes excitability. We found that simulated TTX and RZ blockade of INaP both effectively block intrinsic bursting, either converting bursting dynamics to silence or to tonic spiking, depending in part on the value of gTonic. 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 INaP 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 INaP 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. 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 hNaP 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 Vm- and hNaP-nullclines into the (Vm, hNaP) phase space. The Vm- and hNaP-nullclines are defined from dVm/dt = 0 and dhNaP/dt = 0, respectively. The solution to the latter equation is easy to write down as h n u l l (V m) = h ∞ N a P (V m) = (1 + e h 1 / 2 - V m k) - 1. The Vm-nullcline is determined numerically to form a manifold consisting of three branches that connect at folds or “knees”, forming an N-shape in the (Vm, hNaP) -plane. Intersections between the Vm- and hNaP-nullclines represent equilibrium points of the full system. In the bursting regime, the hNaP-nullcline intersects the middle branch of the Vm-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 Vm- and hNaP-nullclines. Simulated TTX blockade changes the shape of the Vm-nullcline by moving the left and right knees to larger and smaller values of hNaP, respectively, although the right knee is not relevant to the behaviors under study. The former effect corresponds to the increased deinactivation of INaP needed for the neuron to activate after gNaP has been decreased. As gNaP is lowered, the model remains in the bursting regime until the part of the Vm-nullcline near the left knee intersects the hNaP-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 hNaP-nullcline, inducing a shift in the leftward, hyperpolarizing direction. In this case, as hNaP is diminished, bursting continues until the hNaP-nullcline intersects the Vm-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 Fig 3. Although both forms of INaP 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 hNaP rises to a value above that of the left knee of the Vm-nullcline, call it h N a P L K, such that the model’s trajectory can transition to high voltages. The position of the entire Vm-nullcline, including its left knee, depends not only on the extent of INaP blockade but also on the level of input to the neuron. Indeed, hyperpolarizing inputs shift the Vm-nullcline to more hyperpolarized Vm values. As a result, its left branch equilibrium point ends up more hyperpolarized in Vm and, due to the voltage-dependence of INaP, at a more deinactivated, larger hNaP value, call it h N a P E Q. The left knee of the Vm-nullcline also rises to a new value, h ˜ N a P L K > h N a P L K. In theory, as the model settles toward the equilibrium h N a P E Q, although hNaP remains below h ˜ N a P L K, hNaP could exceed the original left knee height, h N a P L K. 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 N a P L K must be less than 1, so that it is physically realizable, (2) the hyperpolarization must be sufficiently strong that h N a P E Q > h N a P L K, 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 INaP in our model. We found that with simulated TTX blockade of INaP rebound bursting is still possible over a narrow range of gNaP, gNaP ∈ (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 N a P L K and thus condition (1), h N a P L K < 1, holds for all levels of Δh1/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–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 INaP, 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 N a P E Q > h N a P L K. If the block is induced with RZ, however, then the trajectory reaches a position in the (Vm, hNaP) 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 INaP even after RZ application. 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 gTonic and the leak reversal potential (ELeak), the latter of which is normally distributed in order to introduce heterogeneity in the population. INaP 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). 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 INaP 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 INaP in respiratory circuits [9,12,13], RZ also inhibits glutamatergic 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), ISynE was blocked in a separate set of simulations by systematically reducing the excitatory synaptic conductance parameter (WMaxE). In contrast to the INaP 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 ∫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 population oscillations before oscillations abruptly stop (Fig 9A). With RZ, however, assuming that both INaP 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 INaP 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 WMaxE decayed exponentially with increasing hyperpolarization of h 1 / 2 N a P (see Materials and methods). To match experimental results of RZ application at 5 μM, WMaxE was maximally reduced by 7% (Fig 9B, Tuning 1), and for RZ at both 10 μM and 20 μM, WMaxE 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 INaP threshold (I N a P T H) required for spike generation and the peak ISynE magnitude during network bursts, and we estimated the hNaP threshold (h N a P T H) 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 ISynE and the firing rate of individual neurons during bursting is a function of the total depolarizing current, determined by both INaP and ISynE. In contrast, network oscillation frequency is determined by the time required for INaP to recover from inactivation to a sufficient threshold (I N a P T H) for spike generation following each network burst. To be more precise, network frequency is a function of the time required for hNaP to recover from inactivation up to a level, h N a P T H, such that I N a P T H is reached following each network burst. For a fixed level of synaptic input, threshold I N a P T H is a constant given by I N a P T H = - g N a P · m N a P · h N a P T H · (V m - E N a) (see Eq 5 in Materials and methods). Thus, as long as synaptic input remains constant, h N a P T H increases as gNaP decreases. In the period of hNaP recovery leading up to burst onset, we indeed have the constant synaptic input ISynE = 0, independent of WMaxE, since the pre-BötC neurons activate synchronously in the regime we are considering. In the case of simulated TTX, while I N a P T H and ISynE did not change, the above reasoning yields an increase in h N a P T H with the progressive decrease in gNaP (Fig 9a1). Note that the increase in the h N a P T H is consistent with Fig 4. Therefore, the reduced network frequency seen with increasing TTX blockade is a direct consequence of an increased hNaP recovery time driven by the increased h N a P T H and the proximity of the nullclines (Fig 4). Because I N a P T H does not vary with gNaP, the level of INaP expressed in bursting neurons is the same as before the blockade. This invariance of INaP together with the fact that TTX does not affect Isyn 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 INaP. In the case of simulated RZ, since gNaP does not vary, h N a P T H 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 hNaP-nullcline), Figs 2 & 4. This shift is accompanied by a reduction of ISynE that follows directly from the reduction of gSynE with simulated RZ. Therefore, changes in network amplitude with simulated RZ arise due to the progressive reduction of ISynE, which results in a decreased firing rate during bursting and de-recruitment of a subset of the neurons. Blockade of INaP 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, INaP 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 INaP 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 INaP 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 INaP 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 INaP block results in a progressive reduction in network amplitude and an increase in frequency before oscillations eventually stop (Fig 10B). The loss of INaP 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 INaP is blocked in pacemaker and non-pacemaker neurons in a random order, then INaP block results in a progressive reduction in network amplitude and frequency before oscillations eventually stop. In all cases, network oscillations stop as soon as INaP is completely blocked in all pacemaker neurons. Finally, to investigate the effects of simulated TTX and RZ block of INaP 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–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 INaP 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 INaP in the intact network, we first characterized the dependence of pre-I population bursting on INaP and synaptic inhibition in the intact network as a function of gTonic (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 INaP = 0 or ISynE = 0. We found that under baseline conditions, bursting in the pre-I population is extremely robust and oscillations continue over all tested values of gTonic, 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 gTonic ∈ (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 ISynE = 0. In contrast, complete blockade of INaP reveals a regime of INaP -dependent bursting for gTonic = 0. 3 − 0. 486, indicated by a transition of the pre-I population from a bursting to a silent mode of activity when INaP = 0. Importantly, these two regimes overlap (for gTonic ∈ (0. 40,0. 486) ), which indicates a region where both INaP 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 INaP during pre-I population bursting, which, as in the pre-I network, modulate INaP inactivation dynamics. As with the pre-I network, there is a threshold level of INaP needed for pre-BötC burst onset, which is realized when hNaP deinactivates to a corresponding threshold level, h N a P T h. One new feature in this case, however, is that h N a P T h depends both on gTonic and on the level of synaptic inhibition to the pre-BötC population, since this parameter and input affect the position of the Vm-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 INaP-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 hNaP, (2) the dynamic range of hNaP over a complete respiratory cycle, and (3) the peak INaP 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 gTonic (Fig 12B, 12C and 12D). Increasing gTonic naturally decreases the hyperpolarization of the pre-I neurons in the expiratory phase (Fig 12B). Due to its effects on the Vm-nullclines for the pre-I neurons, it also lowers I N a P T h and hence h N a P T h. These changes lower the operating range of hNaP (Fig 12C) and correspondingly result in a lower level of INaP 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 N a P T h and hence h N a P T h through the effect of inhibition on the pre-I Vm-nullcline. These changes raise the operating range of hNaP to higher levels (Fig 12C) and correspondingly allow INaP to reach a higher level. Note that the contribution of INaP is small when gTonic is relatively strong and post-I inhibition is weak, whereas the contribution of INaP is large when the relative strengths of gTonic and post-I inhibition are reversed. Interestingly, under the latter conditions, the magnitude of INaP can be larger in the intact network than in the isolated pre-I population due to the post-I population inhibition, which can enhance INaP deinactivation. Finally, we simulated TTX and RZ blockade of INaP in the intact respiratory network (Fig 13). The effects of simulated INaP blockade on the pre-I population dynamics are expected to depend on the value of gTonic 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 gTonic and post-I inhibition for ‘complete’ blockade of INaP by simulated TTX and RZ. Blockade of INaP was considered complete for simulated TTX when gNaP = 0 nS, and for simulated RZ when Δh1/2 = −15 mV (see Fig 2). In the case of simulated TTX, complete blockade of INaP abolishes pre-I population oscillations for gTonic < 0. 486 nS across all levels of post-I inhibition considered (Fig 13A), as expected from the boundary between INaP- and network-dependent bursting regions in Fig 12A. In contrast, for gTonic > 0. 486 nS pre-I population oscillations continue after complete blockade of INaP. In this regime, network amplitude is generally reduced without affecting frequency, except for a frequency reduction for gTonic ≈ 0. 486 nS. The extent of the reduction in amplitude depends on gTonic and the strength of post-I inhibition. The reduction in amplitude is largest when gTonic is weak and post-I inhibition is strong. In the case of simulated RZ, the range of gTonic values where complete blockade of INaP abolishes pre-I population bursting is greatly reduced (Fig 13B) compared to simulated TTX. In general, stronger post-I inhibition decreases the range of gTonic values where INaP block abolishes pre-I population bursting. This relationship is due to the increased recovery of INaP 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 gTonic ≈ 0. 4 nS and post-I inhibition is strong (≳ 0. 5 nS), Fig 13. Importantly, in these simulations, complete blockade of INaP by RZ fails to stop the respiratory rhythm under network conditions where pre-I population bursting is dependent on INaP (Fig 12A). Finally, we compared the effects of simulated RZ blockade of INaP 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 experiments, the maximal 20 μM RZ concentration, which is assumed to completely block INaP, 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 gTonic between 0. 375 nS and 0. 425 nS (Fig 13B and 13C). Interestingly, this range of gTonic values overlaps with the regimes where pre-I population bursting is exclusively INaP -dependent and where bursting is dependent on both INaP and inhibitory network interactions (Fig 12A). Consequently, under these network conditions, our model predicts that selective and complete block of INaP by experimental application of TTX at low concentrations (≤ 20 nM) will abolish the respiratory rhythm in the intact network (Fig 13D).
  This computational study provides novel insight into the role of INaP in respiratory rhythmogenesis and pattern formation. In the intact respiratory network, the role of INaP is not well understood. Current thinking is that in the intact network INaP 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 INaP 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 INaP, however. Moreover, this interpretation assumes that RZ effectively blocks INaP 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 INaP 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 INaP inactivation curve and from simulated TTX blockade based on a reduction in INaP conductance, INaP must strongly contribute to pre-I population bursting. Importantly, this contribution only occurs if, during periods when pre-I neurons are not spiking, INaP strongly deinactivates and hNaP 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 INaP plays a critical role in respiratory rhythmogenesis, which is indicated by the prediction that complete blockade of INaP 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 INaP inactivation. We showed that the dependence of pre-I bursting on INaP 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, INaP 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 INaP 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–20], which, our results show, may compromise the ability of RZ to effectively block INaP. 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 INaP and inhibitory network interactions (Fig 12A & 13). In these simulations, since RZ blockade of INaP fails to stop the rhythm and pre-I bursting is INaP -dependent, this indicates that RZ fails to completely block INaP. In contrast, our simulations predict that under these network conditions a complete blockade of INaP by TTX will stop the respiratory rhythm. Importantly, at low concentrations (≤ 20 nM), TTX has been shown to selectively block INaP 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 INaP in rhythm generation within intact respiratory circuits. One caveat here is that the extent to which rhythm generation in the intact network relies on INaP 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 INaP and a regime where pre-I bursting depends on both INaP 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 INaP -dependent, and complete INaP 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 INaP 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 INaP 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 INaP inhibition. It also did not consider the distinct pharmacological mechanism of action of RZ and INaP block was simulated by reducing gNaP. With the parameter set used for simulations representing the intact network in [13], the pre-I population excitability was set at gtonic = 0. 98 nS and the strength of post-I to pre-I inhibition was 0. 225 nS. In [13] and our simulations, INaP is largely inactivated under these conditions. Consequently, under these parameters, simulated TTX or RZ block of INaP 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, INaP inactivation dynamics, and the role of INaP in rhythm generation in the intact respiratory network may not accurately represent the underlying mechanisms, dynamics, and conditions in these experimental preparations. 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 INaP 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 INaP 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 INaP (Fig 9). The prediction that RZ attenuates ISynE 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 INaP and ISynE in our simulations of the isolated pre-I network, rhythm generation in this system is abolished due to the reduction of INaP, not ISynE. Partial reduction of ISynE alone is not sufficient to stop rhythm generation in our simulations (Fig 8C). These results support the hypothesis that INaP 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 Ca2+ current, inhibition of voltage-gated 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 INaP 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 INaP, 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 off-target drug effects that may arise with bath application. The role of INaP 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 INaP 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 INaP block may provide an explanation for these seemingly contradictory findings. With bath application of TTX or RZ, block of INaP is unlikely to be uniform since drug penetration is dependent on passive diffusion. Therefore, with bath application, INaP 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 INaP 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 Fig 14. Slice thickness varies widely across 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 (non-pacemaker) 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, INaP 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 INaP 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 INaP block. Consistent with this prediction, [25] found pacemaker neurons that are insensitive both to RZ application and to Cd2+ application that blocks Ca2+ 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 INaP 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. In the pre-BötC, coupling of a calcium-activated non-selective cation current (ICAN) and intracellular calcium transients is proposed to represent a biophysical mechanism underlying rhythmogenesis at the cellular- and/or network-level [25,35–41]. For simplicity, however, explicit representation of ICAN and intracellular calcium dynamics were omitted in this study. This choice can be justified by considering the effects of pharmacological block of ICAN in conjunction with recent data-driven computational work. In in vitro slice preparations containing the pre-BötC, pharmacological blockade of ICAN results in a large reduction in the amplitude of network oscillation with no or minor perturbation of frequency, suggesting that the primary role of ICAN in these circuits is amplitude generation/regulation [25,42]. A recent data-driven computational study that included INaP, ICAN, and intracellular calcium dynamics [43] found that in order to match experimental data, ICAN activation must be strongly coupled to synaptically triggered calcium transients, as suggested by [37,39–41,44]. Consequently, ICAN acts as a mechanism that amplifies excitatory synaptic inputs within the pre-BötC. As such, ICAN can be treated as an excitatory post synaptic current that, in our model, could correspond to a portion of the total excitatory current ISynE. Consistent with this idea, blockade of ISynE in our model of the isolated pre-BötC is consistent with the effects seen with experimental blockade of ICAN (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 INaP. 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 INaP -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 INaP block on neuronal activity in these models. For example, inclusion of a Na+/K+ ATPase electrogenic pump in [46], in addition to INaP, dramatically increases the dynamic range where INaP 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–53], but consideration of these additional features is beyond the scope of this study and is left for future work. INaP is not unique to respiratory circuits. Experimental and computational results have suggested an active role for INaP in shaping activity patterns of putative locomotor CPG neurons [54–56] and in contributing to locomotor rhythm generation [55,57,58]. Similar roles are supported for INaP in generating the neural rhythm underlying mastication [59,60]. INaP 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–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 voltage-dependent 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. 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 INaP 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 INaP, and (3) experimental TTX application in the intact network is predicted to terminate respiratory rhythmicity. These findings suggest that INaP 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 INaP 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.
  Neurons were simulated with single compartment models incorporating Hodgkin-Huxley style conductances based on previously described models [13,14,46]. The membrane potential Vm for each neuron is given by the following differential equation: C m d V m d t = I N a + I K + I L e a k + I N a P + I C a + I K C a + I S y n E + I S y n I (1) where Cm = 36. 0 pF is the cell capacitance, INa, IK, ILeak, INaP, ICa, IKCa, ISynE, and ISynI 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 N a = g N a · m N a 3 · h N a · (V m - E N a) (2) I K = g K · m K 4 · (V m - E K) (3) I L e a k = g L e a k · (V m - E L e a k) (4) I N a P = g N a P · m N a P · h N a P · (V m - E N a) (5) I C a = g C a · m C a · h C a · (V m - E C a) (6) I K C a = g K C a · m K C a 2 · (V m - E K) (7) I S y n E = g S y n E · (V m - E S y n E) (8) I S y n I = g S y n I · (V m - E S y n I) (9) where gi is the maximum conductance, Ei is the reversal potential, and mi and hi are gating variables for channel activation and inactivation for current Ii. Values used for gi and Ei are given in Tables 1 and 2. Dynamics of the gating variables m, h for all channels are described by the following differential equation: τ X · d X d t = X ∞ (V) - X; X = { m, h } (10) where X∞ represents steady-state activation/inactivation and τX is a time constant. For INa, INaP, and ICa, the functions X∞ and τX take the forms X ∞ (V) = 1 / (1 + exp (- (V - X 1 / 2) / k X) ), (11) τ X (V) = τ m a x X / cosh ( (V - τ 1 / 2 X) / k τ X). (12) The values of the parameters (X1/2, kX, τ m a x X, τ 1 / 2 X, and k τ X) corresponding to INa INaP and ICa are given in Table 1. For the voltage-gated potassium channel, steady-state activation m ∞ K and time constant τ m K given by the expressions m ∞ K (V) = α ∞ (V) / (α ∞ (V) + β ∞ (V) ), (13) τ m K (V) = 1 / (α ∞ (V) + β ∞ (V) ). (14) where α ∞ (V) = A α · (V + B α) / (1 - exp (- (V + B α) / k α) ), (15) β ∞ (V) = A β · exp (- (V + B β) / k β). (16) The values for the constants Aα, Aβ, Bα, Bβ, kα, and kβ are also given in Table 1. For IKCa, the steady-state activation (m ∞ K C a (C a i n) ) and time constant (τ m K C a (C a i n) ) depend on the intracellular Ca2+ concentration (Cain) as follows: m ∞ K C a (C a i n) = C a i n 2 / (C a i n 2 + β K C a / α K C a), (17) τ m K C a (C a i n) = τ m a x m · 1000 / (α K C a · C a i n 2 + β K C a). (18) The values of αKCa, βKCa and τ m a x m can be found in Table 1. Cain is determined by the balance of Ca2+ influx carried by ICa and efflux via the Ca2+ pump. The dynamics of Cain is described as follows: d C a i n d t = α C a · I C a - (C a i n - C a m i n) τ p u m p (19) where αCa = 2. 5 ⋅ 10−5 mM/fC is a conversion factor relating current to rate of change of Cain, τpump = 500ms is the time constant for the Ca2+ pump and Camin = 5. 0 ⋅ 10−6 mM is a minimal baseline calcium concentration. The total synaptic conductances (g S y n E i and g S y n I i) of the ith target neuron are described by the following equation: g S y n E i = g T o n i c E i + ∑ j ≠ i; n W j, i E · H (t - t j, n) · e - (t - t j, n) / τ s y n E, (20) g S y n I i = ∑ j ≠ i; n W j, i I · H (t - t j, n) · e - (t - t j, n) / τ s y n I. (21) where g T o n i c i is the tonic excitatory conductance, W j, i E and W j, i 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. tj, n is the time at which the nth action potential is generated in neuron j and reaches neuron i. The weights of excitatory and inhibitory conductances were uniformly distributed such that W j, i E = U (0, W M a x E) and W j, i I = U (0, W M a x I) where WMaxE and WMaxI are constants. Parameter values are listed in Tables 1 & 3. 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 Eleak within each population as well as the weights of excitatory (W j, i E) and inhibitory (W j, i I) synaptic connections; see Tables 2 & 3. The tonic excitatory drive (gTonic) 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 (WMaxE) 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), gTonic was treated as a variable to match experimental observations [13]. 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 Vm 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. 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 first-order Euler method with a fixed step-size (Δt) of 0. 025ms.
summary: The application of pharmacological agents that affect transmembrane ionic currents in neurons is a commonly used experimental technique. A simplistic interpretation of experiments involving these agents suggests that antagonist application removes the impacted current and that subsequently observed changes in activity are attributable to the loss of that current’s effects. The more complex reality, however, is that different drugs may have distinct mechanisms of action, some corresponding not to a removal of a current but rather to a changing of its properties. We use computational modeling to explore the implications of the distinct mechanisms associated with two drugs, riluzole and tetrodotoxin, that are often characterized as sodium channel blockers. Through this approach, we offer potential explanations for disparate findings observed in experiments on neural respiratory circuits and show that the experimental results are consistent with a key role for the persistent sodium current in respiratory rhythm generation.
keywords: biotechnology
  knees
  medicine and health sciences
  action potentials
  legs
  neural networks
  engineering and technology
  membrane potential
  electrophysiology
  neuroscience
  pharmacodynamics
  pacemakers
  pharmacology
  body limbs
  neuronal tuning
  bioengineering
  computer and information sciences
  musculoskeletal system
  animal cells
  hyperpolarization
  cellular neuroscience
  cell biology
  anatomy
  physiology
  neurons
  biology and life sciences
  cellular types
  neurophysiology
  medical devices and equipment
title: Effects of persistent sodium current blockade in respiratory circuits depend on the pharmacological mechanism of action and network dynamics