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) 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 inactivation in the hyperpolarizing direction [4-6]. At low concentrations, both TTX (≤ 20 nM) and RZ (≤ 20 μM) have been shown to selectively block I over the fast action-potential generating sodium current (I) [7-10]. A caveat to RZ blockade, however, is that RZ has been shown to inhibit excitatory synaptic transmission at concentrations that affect I [10, 11].

In respiratory circuits, in vitro blocking studies using TTX or RZ have suggested that I 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 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 I block by RZ fails to stop respiratory rhythms in intact preparations [13]. This conclusion is dependent on the assumption that RZ effectively blocks I in an in vivo setting.

This study predicts that after blockade via RZ, I can be reactivated by transient hyperpolarizing perturbations due to the specific mechanism of action of RZ. Simulated TTX blockade of I 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 I 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, while complete block of I 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 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 via RZ and TTX in respiratory neurons of the pre-BötC. Experiments have clearly established the presence of I in neurons of this type and have suggested that neurons exhibiting I-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 -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. Depending on their level of excitability (controlled by g), 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 is required to depolarize neurons above the spiking threshold. During the spiking within a burst, however, I begins to inactivate, which leads to hyperpolarization and burst termination. In the subsequent silent period or interburst interval, the membrane potential V is hyperpolarized, which allows for the voltage-dependent recovery of I from inactivation. After sufficient recovery, sub-threshold activation of I again depolarizes V above the spiking threshold, causing burst initiation and completion of one burst cycle (Fig 1B).

Fig 1

I-dependent intrinsic bursting in a pre-BötC neuron model.

(A) Voltage trace illustrating the dependence of simulated neuronal activity patterns on excitability. Inset (bottom left) shows burst shape. Excitability was varied by linearly increasing the tonic excitatory synaptic conductance (g) from 0.25 nS to 0.45 nS over the course of the simulation (bottom grey ramp). (B) Phase plane plot illustrating the relationship between the slow, voltage-dependent I inactivation variable h and the membrane potential V during the bursting activity displayed in (A). Other model variables are not displayed in this view.

Simulated TTX and RZ blockade of I in steady-state conditions

Bursting in this model is dependent on I, therefore, we first characterized the effects of simulated TTX and RZ blockade on I under steady-state conditions (Fig 2). Since TTX directly obstructs the pores of sodium channels, TTX blockade was simulated by reducing the conductance g of I. Since TTX only affects conductance, the steady-state activation and inactivation parameters of I were not varied (Fig 2A1). We observed that the steady-state current is reduced proportionally to the decrease of g (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 I in steady-state conditions (Fig 2B2). Interestingly, although TTX and RZ blockade of I work through distinct mechanisms, the effects on the peak steady-state current are remarkably similar (Fig 2C).

Fig 2

Effects of simulated TTX and RZ blockade on I under steady-state conditions.

(A1 & B1) Effects of simulated TTX and RZ blockade on voltage-dependent steady-state activation (m∞(V)) and inactivation (h∞(V)) for I. Notice that simulated RZ blockade shifts inactivation in the hyperpolarizing direction and simulated TTX blockade induces no change. (A2 & B2) Effect of simulated TTX and RZ blockade on the I-V curves for I. Notice that simulated TTX and RZ blockade have nearly indistinguishable effects. (C) Peak I current as a function of the extent of simulated TTX and RZ blockade. TTX and RZ blockade are simulated by reducing I conductance (g) and shifting its half-inactivation (h1/2) in the hyperpolarizing direction, respectively.

Simulated TTX and RZ blockade of INaP abolish intrinsic bursting in an isolated model pre-BötC neuron

The activity pattern exhibited by a neuron capable of I-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, represented by reduction of g 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 (g), which tunes excitability. We found that simulated TTX and RZ blockade of I both effectively block intrinsic bursting, either converting bursting dynamics to silence or to tonic spiking, depending in part on the value of g. 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).

Fig 3

Activity patterns in a model pre-BötC neuron as a function of tonic excitatory synaptic drive (g) and (A) level of g and (B) Δh1/2, which are reduced to simulate TTX and RZ block of I, respectively.

Burst frequency (top) and burst duration (bottom) are indicated by color within the bursting regions. Both TTX and RZ effectively abolish intrinsic bursting (colored region): bursting does not occur for g < 2.7nS nor for Δh1/2 < −8.0mV.

Since simulated TTX and RZ blockade have remarkably similar effects on I 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 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 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- and h-nullclines into the (V, h) phase space. The V- and h-nullclines are defined from dV/dt = 0 and dh/dt = 0, respectively. The solution to the latter equation is easy to write down as

The V-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, h)-plane. Intersections between the V- and h-nullclines represent equilibrium points of the full system. In the bursting regime, the h-nullcline intersects the middle branch of the V-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- and h-nullclines. Simulated TTX blockade changes the shape of the V-nullcline by moving the left and right knees to larger and smaller values of h, respectively, although the right knee is not relevant to the behaviors under study. The former effect corresponds to the increased deinactivation of I needed for the neuron to activate after g has been decreased. As g is lowered, the model remains in the bursting regime until the part of the V-nullcline near the left knee intersects the h-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-nullcline, inducing a shift in the leftward, hyperpolarizing direction. In this case, as h is diminished, bursting continues until the h-nullcline intersects the V-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.

Fig 4

Phase plane representation of simulated TTX and RZ block of I in an isolated pre-BötC neuron.

(A) V- and h-nullclines (black and gray, respectively) projected to the (V, h)-plane under control conditions (g = 5.0 nS) in a bursting regime, along with a projection of the bursting trajectory (blue). Note that this view does not represent the I activation variable m. (B) Increasing TTX blockade of I moves the left knee of the V-nullcline to larger values of h, which also abolishes bursting by creating a stable fixed point on the left branch of the V-nullcline via a SNIC (saddle-node on an invariant circle) bifurcation. (C) Increasing RZ blockade of I shifts the h-nullcline to the left, abolishing intrinsic bursting by creating a stable fixed point left branch of the V-nullcline, also via a SNIC bifurcation. g = 0.35 nS in all panels.

Transient hyperpolarization elicits rebound bursting after I blockade

Although both forms of I 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 rises to a value above that of the left knee of the V-nullcline, call it , such that the model’s trajectory can transition to high voltages. The position of the entire V-nullcline, including its left knee, depends not only on the extent of I blockade but also on the level of input to the neuron. Indeed, hyperpolarizing inputs shift the V-nullcline to more hyperpolarized V values. As a result, its left branch equilibrium point ends up more hyperpolarized in V and, due to the voltage-dependence of I, at a more deinactivated, larger h value, call it . The left knee of the V-nullcline also rises to a new value, . In theory, as the model settles toward the equilibrium , although h remains below , h could exceed the original left knee height, . 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) must be less than 1, so that it is physically realizable, (2) the hyperpolarization must be sufficiently strong that , 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 in our model. We found that with simulated TTX blockade of I rebound bursting is still possible over a narrow range of g, g ∈ (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 and thus condition (1), , holds for all levels of Δh1/2.

Fig 5

Hyperpolarized holding potential (V) and duration (color-coded) required to elicit rebound bursting as a function of simulated (A) TTX or (B) RZ blockade of I.

Here, g = 0.35 nS.

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 I, 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), . If the block is induced with RZ, however, then the trajectory reaches a position in the (V, h) 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 even after RZ application.

Fig 6

A transient hyperpolarizing perturbation fails to elicit rebound bursting after simulated (A) TTX but not (B) RZ blockade of I.

(A & B) Top Panels: V- and h-nullclines during baseline, hyperpolarization, and release from hyperpolarization. Arrows indicate the direction of the trajectory. (A & B) Bottom Panels: voltage trace of V during baseline hyperpolarization and release. Initial conditions that generate transient bursting trajectories are indicated by the filled red regions and only appear in the RZ case. For simulated TTX blockade in (A) g = 1.25 nS. In this case, after hyperpolarization (middle panel), the equilibrium point (open circle) lies at an h value that is too low to allow clearance of the left knee after release (right panel). For simulated RZ blockade in (B), Δh1/2 = −12 mV. Here, in the red region, and hence any hyperpolarization that allows the trajectory to enter the red region (e.g., center panel) will yield a burst upon release (right panel). For both simulations, g = 0.35 nS.

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 and the leak reversal potential (E), the latter of which is normally distributed in order to introduce heterogeneity in the population. I 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).

Fig 7

Simulated pre-BötC network activity.

(top) Raster plots showing (A) intrinsic bursting in a subset of neurons in the synaptically uncoupled network and (B) synchronized network bursts in the synaptically coupled network. (bottom) Integrated population activity for the (C) uncoupled and (D) coupled network.

Uniform TTX and RZ block of I 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 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 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), I was blocked in a separate set of simulations by systematically reducing the excitatory synaptic conductance parameter (W). In contrast to the I 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].

Fig 8

Effects of simulated (A) TTX and (B) RZ block of I as well as (C) RZ block of I on network amplitude and frequency of network oscillations in the pre-BötC pre-I neuron population.

Fig 9

Comparison of experimental (colored) and simulated (black) effects of (A) TTX and (B) RZ application in the pre-BötC on the amplitude and frequency of pre-I network oscillations.

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 and excitatory synapses whereas TTX only affects I. In the simulations, the relevant values at the end points (where network oscillations stop) are as follows: (A) g = 3.52; (B) (top trace) Δh1/2 = −5.0 mV, W = 0.0255 nS, (bottom trace) Δh1/2 = −4.5 mV, W = 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, peak I, and the I inactivation threshold () required to initiate bursting. was defined as the maximal value of the mean population h prior to burst initiation.

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 I 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 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 decayed exponentially with increasing hyperpolarization of (see Materials and methods). To match experimental results of RZ application at 5 μM, W was maximally reduced by 7% (Fig 9B, Tuning 1), and for RZ at both 10 μM and 20 μM, W 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 threshold () required for spike generation and the peak I magnitude during network bursts, and we estimated the h threshold () 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 and the firing rate of individual neurons during bursting is a function of the total depolarizing current, determined by both I and I. In contrast, network oscillation frequency is determined by the time required for I to recover from inactivation to a sufficient threshold () for spike generation following each network burst. To be more precise, network frequency is a function of the time required for h to recover from inactivation up to a level, , such that is reached following each network burst. For a fixed level of synaptic input, threshold is a constant given by (see Eq 5 in Materials and methods). Thus, as long as synaptic input remains constant, increases as g decreases. In the period of h recovery leading up to burst onset, we indeed have the constant synaptic input I = 0, independent of W, since the pre-BötC neurons activate synchronously in the regime we are considering.

In the case of simulated TTX, while and I did not change, the above reasoning yields an increase in with the progressive decrease in g (Fig 9a1). Note that the increase in the is consistent with Fig 4. Therefore, the reduced network frequency seen with increasing TTX blockade is a direct consequence of an increased h recovery time driven by the increased and the proximity of the nullclines (Fig 4). Because does not vary with g, the level of I expressed in bursting neurons is the same as before the blockade. This invariance of I together with the fact that TTX does not affect I 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.

In the case of simulated RZ, since g does not vary, 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-nullcline), Figs 2 & 4. This shift is accompanied by a reduction of I that follows directly from the reduction of g with simulated RZ. Therefore, changes in network amplitude with simulated RZ arise due to the progressive reduction of I, 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 in the pre-BötC pre-I network

Blockade of I 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 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 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 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 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 block results in a progressive reduction in network amplitude and an increase in frequency before oscillations eventually stop (Fig 10B). The loss of I 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 is blocked in pacemaker and non-pacemaker neurons in a random order, then I 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 is completely blocked in all pacemaker neurons.

Fig 10

Simulated progressive nonuniform I block in the isolated pre-I network.

(A,B,C) From left to right: network burst with neurons color-coded as pacemakers (red) or non-pacemakers (blue) and numbered based on order affected from 1 (first affected) to 50 (last affected); amplitude and frequency of isolated pre-I population oscillations as a function of the percentage of the pre-I population where I is completely blocked. For a given neuron, I is considered completely blocked when g = 0 nS in the case of TTX (A1,B1,C1) or when Δh1/2 = −15 mV and I is attenuated by 25% for RZ (A2,B2,C2), see Figs 2 and 8. Error bars in C1 and C2 indicate the ±SEM of ten trials where I is progressively blocked across the network in random order.

Simulated TTX and RZ block of I in the intact respiratory network

Finally, to investigate the effects of simulated TTX and RZ block of I 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.

Fig 11

Simulated intact respiratory network.

(A) Circuit diagram of the intact respiratory network composed of the inhibitory (post-I, aug-E early-I) and excitatory (pre-I) subpopulations. (B) Integrated subpopulation spiking activity. Amplitudes of each subpopulation were normalized. (C) Spiking in example neurons from each subpopulation.

In the intact respiratory network, the respective contributions of I 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 in the intact network, we first characterized the dependence of pre-I population bursting on I and synaptic inhibition in the intact network as a function of g (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 = 0 or I = 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, 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 ∈ (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 = 0. In contrast, complete blockade of I reveals a regime of I -dependent bursting for g = 0.3 − 0.486, indicated by a transition of the pre-I population from a bursting to a silent mode of activity when I = 0. Importantly, these two regimes overlap (for g ∈ (0.40, 0.486)), which indicates a region where both I and inhibitory synaptic interactions are critical for rhythm generation/bursting in the pre-I population (Fig 12A).

Fig 12

Characterization of the intact respiratory network behavior.

(A) Identification of I and/or network dependent regimes of pre-I population bursting as a function of tonic excitatory drive (g). Effect of post-I inhibition strength and g on (B) post-inspiratory phase hyperpolarization, (C) h dynamic range, and (D) the peak inspiratory phase I in the pre-I population. Notice that the I level is strongly affected by the magnitude of post-I inhibition and g. Strong post-I inhibition increases the magnitude of post-inspiratory phase hyperpolarization, which decreases I inactivation and increases the peak I. In contrast, strong g decreases post-inspiratory phase hyperpolarization, which increases I inactivation and decreases the peak I. Dashed lines in C indicate the dynamic range of h in the isolated pre-I population under baseline conditions used in Fig 7.

Next we characterized network parameters affecting the magnitude/contribution of I during pre-I population bursting, which, as in the pre-I network, modulate I inactivation dynamics. As with the pre-I network, there is a threshold level of I needed for pre-BötC burst onset, which is realized when h deinactivates to a corresponding threshold level, . One new feature in this case, however, is that depends both on g and on the level of synaptic inhibition to the pre-BötC population, since this parameter and input affect the position of the V-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-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, (2) the dynamic range of h over a complete respiratory cycle, and (3) the peak I 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 (Fig 12B, 12C and 12D).

Increasing g naturally decreases the hyperpolarization of the pre-I neurons in the expiratory phase (Fig 12B). Due to its effects on the V-nullclines for the pre-I neurons, it also lowers and hence . These changes lower the operating range of h (Fig 12C) and correspondingly result in a lower level of I 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 and hence through the effect of inhibition on the pre-I V-nullcline. These changes raise the operating range of h to higher levels (Fig 12C) and correspondingly allow I to reach a higher level. Note that the contribution of I is small when g is relatively strong and post-I inhibition is weak, whereas the contribution of I is large when the relative strengths of g and post-I inhibition are reversed. Interestingly, under the latter conditions, the magnitude of I 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 deinactivation.

Finally, we simulated TTX and RZ blockade of I in the intact respiratory network (Fig 13). The effects of simulated I blockade on the pre-I population dynamics are expected to depend on the value of g 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 and post-I inhibition for ‘complete’ blockade of I by simulated TTX and RZ. Blockade of I was considered complete for simulated TTX when g = 0 nS, and for simulated RZ when Δh1/2 = −15 mV (see Fig 2).

Fig 13

Predicted effect of blockade of I in the intact network.

Relative effects of “complete” blockade of I by simulated (A) TTX and (B) RZ application on the amplitude (left) and frequency (right) as a function of tonic excitatory drive (g) and the strength of post-I inhibition. Oval shape indicates parameters capable of matching the maximal relative changes in amplitude and frequency seen with experimental application of RZ. The points labeled C and D indicate the g and post-I inhibition values used in panels C and D. (C) Comparison of experimental and simulated RZ blockade of I. Simulated RZ application closely matches experimental data when g = 0.4 and post-I inhibition = 0.82. Experimental data is adapted from [13]. (D) Predicted effect of experimental TTX blockade of I in the intact network under identical conditions to panel C.

In the case of simulated TTX, complete blockade of I abolishes pre-I population oscillations for g < 0.486 nS across all levels of post-I inhibition considered (Fig 13A), as expected from the boundary between I- and network-dependent bursting regions in Fig 12A. In contrast, for g > 0.486 nS pre-I population oscillations continue after complete blockade of I. In this regime, network amplitude is generally reduced without affecting frequency, except for a frequency reduction for g ≈ 0.486 nS. The extent of the reduction in amplitude depends on g and the strength of post-I inhibition. The reduction in amplitude is largest when g is weak and post-I inhibition is strong.

In the case of simulated RZ, the range of g values where complete blockade of I 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 values where I block abolishes pre-I population bursting. This relationship is due to the increased recovery of I 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 ≈ 0.4 nS and post-I inhibition is strong (≳ 0.5 nS), Fig 13. Importantly, in these simulations, complete blockade of I by RZ fails to stop the respiratory rhythm under network conditions where pre-I population bursting is dependent on I (Fig 12A).

Finally, we compared the effects of simulated RZ blockade of I 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 I, 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 between 0.375 nS and 0.425 nS (Fig 13B and 13C). Interestingly, this range of g values overlaps with the regimes where pre-I population bursting is exclusively I -dependent and where bursting is dependent on both I and inhibitory network interactions (Fig 12A). Consequently, under these network conditions, our model predicts that selective and complete block of I 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 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 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 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 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 blockade in the intact respiratory network predict that experimental application of RZ, but not TTX, will fail to effectively block I and respiratory rhythmicity.

Insights into the role of I in the intact respiratory network

This computational study provides novel insight into the role of I in respiratory rhythmogenesis and pattern formation. In the intact respiratory network, the role of I is not well understood. Current thinking is that in the intact network I 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 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, however. Moreover, this interpretation assumes that RZ effectively blocks I 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 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 inactivation curve and from simulated TTX blockade based on a reduction in I conductance, I must strongly contribute to pre-I population bursting. Importantly, this contribution only occurs if, during periods when pre-I neurons are not spiking, I strongly deinactivates and h 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 plays a critical role in respiratory rhythmogenesis, which is indicated by the prediction that complete blockade of I 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 inactivation.

We showed that the dependence of pre-I bursting on I 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 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 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 I. 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 and inhibitory network interactions (Fig 12A & 13). In these simulations, since RZ blockade of I fails to stop the rhythm and pre-I bursting is I -dependent, this indicates that RZ fails to completely block I. In contrast, our simulations predict that under these network conditions a complete blockade of I by TTX will stop the respiratory rhythm. Importantly, at low concentrations (≤ 20 nM), TTX has been shown to selectively block I 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 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 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 and a regime where pre-I bursting depends on both I 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 -dependent, and complete I 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 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 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 inhibition. It also did not consider the distinct pharmacological mechanism of action of RZ and I block was simulated by reducing g. With the parameter set used for simulations representing the intact network in [13], the pre-I population excitability was set at g = 0.98 nS and the strength of post-I to pre-I inhibition was 0.225 nS. In [13] and our simulations, I is largely inactivated under these conditions. Consequently, under these parameters, simulated TTX or RZ block of I 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 inactivation dynamics, and the role of I 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 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 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 (Fig 9). The prediction that RZ attenuates I 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 and I in our simulations of the isolated pre-I network, rhythm generation in this system is abolished due to the reduction of I, not I. Partial reduction of I alone is not sufficient to stop rhythm generation in our simulations (Fig 8C). These results support the hypothesis that I 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 I 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, 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.

Non-uniform pharmacological blockade in the pre-BötC pre-I population

The role of I 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 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 block may provide an explanation for these seemingly contradictory findings.

With bath application of TTX or RZ, block of I is unlikely to be uniform since drug penetration is dependent on passive diffusion. Therefore, with bath application, I 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 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.

Fig 14

Sketch of proposed distribution of pacemaker and non-pacemaker neurons within the pre-BötC as well as the the direction and depth of TTX and RZ penetration in thick and thin in vitro slice preparation containing the pre-BötC.

The direction of diffusion and maximal diffusion depth in the thin and thick slice sketches are indicated by arrows and vertical dashed lines. Notice that in a thick slice non-pacemaker neurons are affected first and TTX and RZ fail to fully penetrate the thick slice.

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 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 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 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 I 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) 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 I and intracellular calcium dynamics were omitted in this study. This choice can be justified by considering the effects of pharmacological block of I in conjunction with recent data-driven computational work. In in vitro slice preparations containing the pre-BötC, pharmacological blockade of I 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 in these circuits is amplitude generation/regulation [25, 42]. A recent data-driven computational study that included I, I, and intracellular calcium dynamics [43] found that in order to match experimental data, I activation must be strongly coupled to synaptically triggered calcium transients, as suggested by [37, 39–41, 44]. Consequently, I acts as a mechanism that amplifies excitatory synaptic inputs within the pre-BötC. As such, I can be treated as an excitatory post synaptic current that, in our model, could correspond to a portion of the total excitatory current I. Consistent with this idea, blockade of I in our model of the isolated pre-BötC is consistent with the effects seen with experimental blockade of I (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. 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 -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 block on neuronal activity in these models. For example, inclusion of a Na+/K+ ATPase electrogenic pump in [46], in addition to I, dramatically increases the dynamic range where I 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.

Broader implications

I is not unique to respiratory circuits. Experimental and computational results have suggested an active role for I 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 I in generating the neural rhythm underlying mastication [59, 60]. I 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.

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 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, and (3) experimental TTX application in the intact network is predicted to terminate respiratory rhythmicity. These findings suggest that I 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 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.

Materials and methods

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 for each neuron is given by the following differential equation: where C = 36.0 pF is the cell capacitance, I, I, I, I, I, I, I, and I 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: where g is the maximum conductance, E is the reversal potential, and m and h are gating variables for channel activation and inactivation for current I. Values used for g and E are given in Tables 1 and 2.

Table 1

Ionic channel parameters.

U(a, b) indicates a uniform distribution from a to b.

Channel	Parameters
INa	ENa = 55.0 mV
	m1/2 = −43.8 mV	km = 6.0 mV	τmaxm=0.25ms	τ1/2m=-43.8mV	kτm=14.0mV
	h1/2 = −67.5 mV	kh = −10.8 mV	τmaxh=8.46ms	τ1/2h=-67.5mV	kτh=12.8mV
IK	EK = −94.4 mV
	Aα = 0.01	Bα = 44.0 mV	kα = 5.0 mV
	Aβ = 0.17	Bβ = 49.0 mV	kβ = 40.0 mV
ILeak	ELeak = −68 mV (Isolated Pre-I neuron)
INaP	m1/2 = −47.1 mV	km = 3.1 mV	τmaxm=1.0ms	τ1/2m=-47.1mV	kτm=6.2mV
	h1/2 = −60.0 mV	kh = −9.0 mV	τmaxh=5000ms	τ1/2h=-60.0mV	kτh=9.0mV
ICa	ECa = 13.27 ⋅ ln(Caout/Cain)		Caout = 4.0 mM
	m1/2 = −27.5 mV	km = 5.7 mV	τm = 0.5 ms
	h1/2 = −52.4 mV	kh = −5.2 mV	τh = 18.0 ms
IKCa	αKCa = 1.25 ⋅ 108	βKCa = 2.5	τmaxm=U(1,10)ms
ISynE	ESynE = 0.0 mV	τSynE = 5.0 ms
ISynI	ESynI = −75 mV	τSynI = 15.0 ms

Table 2

Neuronal type specific parameters.

U(a, b) indicates a uniform distribution from a to b.

Type	gNa (nS)	gK (nS)	gLeak (nS)	ELeak (mV)	gNaP (nS)	gCa (nS)	gKCa (nS)	gTonicE (nS)
Pre-I	170	180	2.25	U(−66.5, −69.5)	5.0	-	-	0 − 0.5
Early-I	400	250	6.0	U(−58.8, −61.2)	-	0.34	U(3.0, 6.0)	0.25
Aug-E	400	250	6.0	U(−58.8, −61.2)	-	0.042	U(3.0, 6.0)	0.475
Post-I	400	250	6.0	U(−58.8, −61.2)	-	0.17	U(3.0, 6.0)	0.3

Dynamics of the gating variables m, h for all channels are described by the following differential equation: where X∞ represents steady-state activation/inactivation and τ is a time constant. For I, I, and I, the functions X∞ and τ take the forms The values of the parameters (X1/2, k, , , and ) corresponding to I I and I are given in Table 1.

For the voltage-gated potassium channel, steady-state activation and time constant given by the expressions where The values for the constants A, A, B, B, k, and k are also given in Table 1.

For I, the steady-state activation () and time constant () depend on the intracellular Ca2+ concentration (Ca) as follows: The values of α, β and can be found in Table 1.

Ca is determined by the balance of Ca2+ influx carried by I and efflux via the Ca2+ pump. The dynamics of Ca is described as follows: where α = 2.5 ⋅ 10−5 mM/fC is a conversion factor relating current to rate of change of Ca, τ = 500ms is the time constant for the Ca2+ pump and Ca = 5.0 ⋅ 10−6 mM is a minimal baseline calcium concentration.

The total synaptic conductances ( and ) of the i target neuron are described by the following equation: where is the tonic excitatory conductance, and 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. τ and τ are exponential decay constants for excitatory and inhibitory synapses. t is the time at which the n action potential is generated in neuron j and reaches neuron i. The weights of excitatory and inhibitory conductances were uniformly distributed such that and where W and W are constants. Parameter values are listed in Tables 1 & 3.

Table 3

Maximal weight of excitatory (W) and inhibitory (W) synaptic connections in nS.

		WMaxE Target (i)				WMaxI Target (i)
		Pre-I	Early-I	Aug-E	Post-I	Pre-I	Early-I	Aug-E	Post-I
source (j)	Pre-I	0.03	0.125	-	-	-	-	-	-
	Early-I	-	-	-	-	-	-	0.22	0.02635
	Aug-E	-	-	-	-	0.0075	0.05	-	0.01125
	Post-I	-	-	-	-	0.75	0.372	0.15	-

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 within each population as well as the weights of excitatory () and inhibitory () synaptic connections; see Tables 2 & 3. The tonic excitatory drive (g) 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) 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 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 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 first-order Euler method with a fixed step-size (Δt) of 0.025ms.