1. Computational and Systems Biology
  2. Neuroscience
Download icon

Biophysical mechanisms in the mammalian respiratory oscillator re-examined with a new data-driven computational model

  1. Ryan S Phillips
  2. Tibin T John
  3. Hidehiko Koizumi
  4. Yaroslav I Molkov
  5. Jeffrey C Smith  Is a corresponding author
  1. National Institute of Neurological Disorders and Stroke, National Institutes of Health, United States
  2. University of New Hampshire, United States
  3. Georgia State University, United States
Research Article
  • Cited 1
  • Views 477
  • Annotations
Cite this article as: eLife 2019;8:e41555 doi: 10.7554/eLife.41555

Abstract

An autorhythmic population of excitatory neurons in the brainstem pre-Bötzinger complex is a critical component of the mammalian respiratory oscillator. Two intrinsic neuronal biophysical mechanisms—a persistent sodium current (INaP) and a calcium-activated non-selective cationic current (ICAN)—were proposed to individually or in combination generate cellular- and circuit-level oscillations, but their roles are debated without resolution. We re-examined these roles in a model of a synaptically connected population of excitatory neurons with ICAN and INaP. This model robustly reproduces experimental data showing that rhythm generation can be independent of ICAN activation, which determines population activity amplitude. This occurs when ICAN is primarily activated by neuronal calcium fluxes driven by synaptic mechanisms. Rhythm depends critically on INaP in a subpopulation forming the rhythmogenic kernel. The model explains how the rhythm and amplitude of respiratory oscillations involve distinct biophysical mechanisms.

https://doi.org/10.7554/eLife.41555.001

Introduction

Defining cellular and circuit mechanisms generating the vital rhythm of breathing in mammals remains a fundamental unsolved problem of wide-spread interest in neurophysiology (Richter and Smith, 2014; Del Negro et al., 2018; Ramirez and Baertsch, 2018a), with potentially far-reaching implications for understanding mechanisms of oscillatory circuit activity and rhythmic motor pattern generation in neural systems (Marder and Calabrese, 1996; Buzsaki, 2006; Grillner, 2006; Kiehn, 2006). The brainstem pre-Bötzinger complex (pre-BötC) region (Smith et al., 1991) located in the ventrolateral medulla oblongata is established to contain circuits essential for respiratory rhythm generation (Smith et al., 2013; Del Negro et al., 2018), but the operational cellular biophysical and circuit synaptic mechanisms are continuously debated. Pre-BötC excitatory neurons and circuits have autorhythmic properties and drive motor circuits that can be isolated and remain rhythmically active in living rodent brainstem slices in vitro. Numerous experimental and theoretical analyses have focused on the rhythmogenic mechanisms operating in these in vitro conditions to provide insight into biophysical and circuit processes involved, with potential relevance for rhythm generation during breathing in vivo (Feldman and Del Negro, 2006; Lindsey et al., 2012; Richter and Smith, 2014; Ramirez and Baertsch, 2018b). The ongoing rhythmic activity in vitro has been suggested to arise from a variety of cellular and circuit biophysical mechanisms including from a subset(s) of intrinsically bursting neurons which, through excitatory synaptic interactions, recruit and synchronize neurons within the network (pacemaker-network models) (Butera et al., 1999b; Ramirez et al., 2004; Toporikova and Butera, 2011; Chevalier et al., 2016), or as an emergent network property through recurrent excitation (e.g. Rekling and Feldman, 1998; Jasinski et al., 2013) and/or synaptic depression (group pacemaker model) (Rubin et al., 2009a; Del Negro et al., 2010).

From these previous analyses, involvement of two possible cellular-level biophysical mechanisms have been proposed. One based on a slowly inactivating persistent sodium current (INaP) (Butera et al., 1999a), and the other on a calcium-activated non-selective cation current (ICAN) coupled to intracellular calcium (Cai) dynamics (for reviews see Rybak et al., 2014; Del Negro et al., 2010), or a combination of both mechanisms (Thoby-Brisson and Ramirez, 2001Jasinski et al., 2013; Peña et al., 2004). Despite the extensive experimental and theoretical investigations of these sodium- and calcium-based mechanisms, the actual roles of INaP, ICAN and the critical source(s) of Cai transients in the pre-BötC are still unresolved. Furthermore, in pre-BötC circuits, the process of rhythm generation must be associated with an amplitude of circuit activity sufficient to drive downstream circuits to produce adequate inspiratory motor output. Biophysical mechanisms involved in generating the amplitude of pre-BötC circuit activity have also not been established.

INaP is proposed to mediate an essential oscillatory burst-generating mechanism since pharmacologically inhibiting INaP abolishes intrinsic neuronal rhythmic bursting as well as pre-BötC circuit inspiratory activity and rhythmic inspiratory motor output in vitro (Koizumi et al., 2008Toporikova et al., 2015), although some studies suggest that block of both INaP and ICAN are necessary to disrupt rhythmogenesis in vitro (Peña et al., 2004). Theoretical models of cellular and circuit activity based on INaP-dependent bursting mechanisms closely reproduce experimental observations such as voltage-dependent frequency control, spike-frequency adaptation during bursts, and pattern formation of inspiratory motor output (Butera et al., 1999b; Pierrefiche et al., 2004; Smith et al., 2007). This indicates the plausibility of INaP-dependent rhythm generation.

In the pre-BötC, ICAN was originally postulated to underlie intrinsic pacemaker-like oscillatory bursting at the cellular level and contribute to circuit-level rhythm generation, since intrinsic bursting in a subset of neurons in vitro was found to be terminated by the ICAN inhibitor flufenamic acid (FFA) (Peña et al., 2004). Furthermore, inhibition of ICAN in the pre-BötC reduces the amplitude of the rhythmic depolarization (inspiratory drive potential) driving neuronal bursting and can eliminate inspiratory motor activity in vitro (Pace et al., 2007). ICAN became the centerpiece of the 'group pacemaker' model for rhythm generation, in which this conductance was proposed to be activated by inositol trisphosphate (IP3) receptor/ER-mediated intracellular calcium fluxes initiated via glutamatergic metabotropic receptor-mediated signaling in the pre-BötC excitatory circuits (Del Negro et al., 2010). The molecular correlate of ICAN was postulated to be the transient receptor potential channel M4 (TRPM4) (Mironov, 2008; Pace et al., 2007), one of the two known Ca2+-activated TRP channels (Guinamard et al., 2010; Ullrich et al., 2005), or alternatively, by the transient receptor potential channels C3/7 (TRPC3/7) (Ben-Mabrouk and Tryba, 2010); however, these latter channels are not known to be Ca2+-activated (Clapham, 2003). TRPM4 and TRPC3 have now been identified by immunolabeling and RNA expression profiling in pre-BötC inspiratory neurons in vitro (Koizumi et al., 2018).

Investigations into the sources of intracellular Ca2+ activating ICAN/TRPM4 suggested that (1) somatic calcium transients from voltage-gated sources do not contribute to the inspiratory drive potential (Morgado-Valle et al., 2008), (2) IP3/ER-mediated intracellular Ca2+ release does not contribute to inspiratory rhythm generation in vitro, and (3) in the dendrites calcium transients may be triggered by excitatory synaptic inputs and travel in a wave propagated to the soma (Mironov, 2008). Theoretical studies have demonstrated the plausibility of Cai-ICAN-dependent bursting (Rubin et al., 2009b; Toporikova and Butera, 2011); however, these models omit INaP and/or depend on additional unproven mechanisms to generate intracellular calcium oscillations to provide burst termination, such as IP3-dependent calcium-induced calcium release (Toporikova and Butera, 2011), partial depolarization block of action potentials (Rubin et al., 2009a), and the Na+/K+ pump (Jasinski et al., 2013). Interestingly, pharmacological inhibition of ICAN/TRPM4 has been shown to produce large reductions in the amplitude of pre-BötC inspiratory neuron population activity with essentially no, or minor perturbations of inspiratory rhythm (Peña et al., 2004). These observations constrain the role of ICAN, and require theoretical re-examination of pre-BötC neuronal conductance mechanisms and network dynamics, particularly how rhythm generation mechanisms can be independent of ICAN-dependent mechanisms that regulate the amplitude of network activity.

In this theoretical study, we examine the role of ICAN in pre-BötC excitatory circuits by considering two plausible mechanisms of intracellular calcium fluxes: (1) from voltage-gated and (2) from synaptically activated sources. We deduce that ICAN is primarily activated by calcium transients that are coupled to rhythmic excitatory synaptic inputs originating from INaP-dependent bursting inspiratory neurons. Additionally, we show that ICAN contributes to the inspiratory drive potential by mirroring the excitatory synaptic current. This concept is consistent with a mechanism underlying generation of the inspiratory drive potential involving a synaptic-based ICAN activation described in previous work (Pace et al., 2007; Rubin et al., 2009a). Our model explains the experimental observations obtained from in vitro neonatal rodent slices isolating the pre-BötC, showing large reductions in network activity amplitude by inhibiting ICAN/TRPM4 without perturbations of inspiratory rhythm generation in pre-BötC excitatory circuits in vitro. The model supports the concept that ICAN activation in a subpopulation of pre-BötC excitatory neurons is critically involved in amplifying synaptic drive from a subset of neurons whose rhythmic bursting is critically dependent on INaP and forms the kernel for rhythm generation in vitro. The model suggests how the functions of generating the rhythm and amplitude of inspiratory oscillations in pre-BötC excitatory circuits are determined by distinct biophysical mechanisms.

Results

g-CAN variation has opposite effects on amplitude and frequency of network bursting in the CaV and CaSyn models

Experimental work (Peña et al., 2004) has demonstrated that pharmacological inhibition of ICAN/TRPM4 in the pre-BötC from in vitro neonatal mouse/rat slice preparations, strongly reduces the amplitude of (or completely eliminates) the inspiratory hypoglossal (XII) motor output, as well as the amplitude of pre-BötC excitatory circuit activity that is highly correlated with the decline of XII activity, while having relatively little effect on inspiratory burst frequency. Here, we systematically examine in our model the relationship between ICAN conductance (g-CAN) and the amplitude and frequency of network activity for voltage-gated (CaV) and synaptically activated sources (CaSyn) of intracellular calcium in a heterogeneous network of 100 synaptically coupled single-compartment pre-BötC model excitatory neurons. In addition to voltage-gated and synaptically activated calcium currents, each model neuron incorporates voltage-gated action potential generating currents, as well as ICAN, INaP, leak, and excitatory synaptic currents adapted from the conductance-based biophysical model of Jasinski et al. (2013) (see Materials and methods for a full model description). We note that in our full model, ICAN is activated by both voltage-gated and synaptic mechanisms, consistent with experimental results (Thoby-Brisson and Ramirez, 2001; Pace et al., 2007; Morgado-Valle et al., 2008). Initial separate consideration of the CaV and CaSyn provides a means to deduce the relative contribution of these two general sources of intracellular calcium to ICAN activation. We found that reduction of g-CAN drives opposing effects on network activity amplitude (spike/s/neuron) and frequency that are dependent on the source of intracellular calcium transients (Figure 1). The network activity amplitude is a measure of the average neuronal population firing rate and is defined by the number of spikes generated by the network per 50 ms bin divided by the number of neurons in the network. In the CaV network, where calcium influx is generated exclusively from voltage-gated calcium channels, increasing g-CAN has no effect on amplitude but increases the frequency of network oscillations (Figure 1A, C, D). Conversely, in the CaSyn network where calcium influx is generated exclusively by excitatory synaptic input, increasing g-CAN strongly increases the amplitude and slightly decreases the oscillation frequency (Figure 1B, C, D).

Effects of subthreshold activation of ICAN on network frequency

In INaP-dependent bursting neurons in the pre-BötC, bursting frequency depends on their excitability (i.e. baseline membrane potential) which can be controlled in different ways, for example by directly injecting a depolarizing current (Smith et al., 1991Del Negro et al., 2005; Yamanishi et al., 2018) or varying the conductance and/or reversal potentials of some ionic channels (Butera et al., 1999a). Due to their relatively short duty cycle, the bursting frequency in these neurons is largely determined by the interburst interval, defined as the time between the end of one burst and the start of the next. During the burst, INaP slowly inactivates (Butera et al., 1999a) resulting in burst termination and abrupt neuronal hyperpolarization. The interburst interval is then determined by the amount of time required for INaP to recover from inactivation and return the membrane potential back to the threshold for burst initiation. This process is governed by the kinetics of INaP inactivation gating variable hNaP. Higher neuronal excitability reduces the value of hNaP required to initiate bursting. Consequently, the time required to reach this value is decreased, which results in a shorter interburst interval and increased frequency.

To understand how changing g-CAN affects network bursting frequency, we quantified the values of hNaP averaged over all rhythm-generating pacemaker neurons immediately preceding each network burst and, also, the average ICAN values between the bursts in the CaV and Casyn networks (Figure 2). In the CaV network, ICa as modeled remains residually activated between the bursts thus creating the background calcium concentration which partially activates ICAN. Therefore, between the bursts ICAN functions as a depolarizing leak current. Consistently, we found that in the CaV network increasing g-CAN increases ICAN (Figure 2A), progressively depolarizing the network, which reduces the hNaP threshold for burst initiation (Figure 2B) and, thus, increases network oscillation frequency (Figure 1D).

Manipulations of g-CAN in the CaV and CaSyn networks produce opposite effects on network activity amplitude (spikes/s) and frequency.

(A and B) Histograms of neuronal population activity amplitude in the CaV, and CaSyn networks with linearly increasing g-CAN. (C) Plot of g-CAN (% of the baseline mean value for the simulated population) vs. network activity amplitude for the CaV and CaSyn networks in A and B. (D) Plot of g-CAN (%) vs. network frequency for the CaV and CaSyn networks in A and B. CaV network parameters: g-Ca=1.0 (nS), PCa=0.0, PSyn=0.05 and Wmax=0.2 (nS). CaSyn network parameters: g-Ca=0 (nS), PCa=0.01, PSyn=0.05 and Wmax=0.2 (nS).

https://doi.org/10.7554/eLife.41555.002
Calcium source and g-CAN-dependent effects on cellular properties regulating network frequency for the simulations presented in Figure 1.

(A) Average magnitude of ICAN in pacemaker neurons during the interburst interval for the CaV (red) and CaSyn (blue) networks. (B) Average inactivation (hNaP ) of the burst generating current INaP in pacemaker neurons immediately preceding each network burst as a function of g-CAN (%) for the voltage-gated and synaptic calcium networks. CaV network parameters: g-Ca=1.0 (nS), PCa=0.0, PSyn=0.05 and Wmax=0.2 (nS). CaSyn network parameters: g-Ca=0 (nS), PCa=0.01, PSyn=0.05 and Wmax=0.2 (nS).

https://doi.org/10.7554/eLife.41555.003

In the Casyn model, the intracellular calcium depletes entirely during the interburst interval. Consequently, increasing g-CAN has no effect on ICAN (Figure 2A) and frequency is essentially unaffected (Figure 1D).

Changes in network activity amplitude are driven by recruitment of neurons

As previously stated, the network activity amplitude is defined as the total number of spikes produced by the network per a time bin divided by the number of neurons in the network. Consequently, changes in amplitude can only occur by increasing the number of neurons participating in bursts (recruitment) and/or increasing the firing rate of the recruited neurons. To analyze changes in amplitude, we quantified the number of recruited neurons (Figure 3A) and the average spike frequency in recruited neurons (Figure 3B) as a function of g-CAN for both network models. In the CaV network, increasing g-CAN increases the number of recruited neurons (Figure 3A), but decreases the average spiking frequency in recruited neurons (Figure 3B) which, together result in no change in amplitude (Figure 1C). In the CaSyn network, increasing g-CAN strongly increases the number of recruited neurons (Figure 3A) and increases the spike frequency of recruited neurons (Figure 3B) resulting in a large increase in network activity amplitude (Figure 3C).

Calcium source and g-CAN-dependent effects on cellular properties regulating network activity amplitude for the simulations presented in Figure 1.

(A) Number of recruited neurons in the modeled population of 100 neurons as a function of g-CAN (%) for voltage-gated and synaptic calcium sources. The number of recruited neurons is defined as the peak number of spiking neurons per bin during a network burst. (B) Average spiking frequency of recruited neurons as a function of g-CAN for the voltage-gated and synaptic calcium mechanism. Average spiking frequency is defined the number of spikes per bin divided by the number of recruited neurons. The parameters used in these simulations are: CaVg-Ca=1.0 (nS), PCa=0.0, PSyn=0.05 and Wmax=0.2 (nS)CaSyng-Ca=0 (nS), PCa=0.01, PSyn=0.05 and Wmax=0.2 (nS).

https://doi.org/10.7554/eLife.41555.004

Manipulating g-CAN in the CaSyn model is qualitatively equivalent to changing the strength of synaptic interactions

Since changes in g-CAN in the CaSyn model primarily affect network activity amplitude through recruitment of neurons, and the network activity amplitude strongly depends on the strength of synaptic interactions, we next examined the relationship between g-CAN, synaptic strength and network activity amplitude and frequency (Figure 4). Synaptic strength is defined as the number of neurons multiplied by the synaptic connection probability multiplied by the average weight of synaptic connections (N·PSyn·12Wmax), where the weight of synaptic connections Wmax ranges from 0.0 to 1.0 nS. We found that the effects of varying g-CAN or the synaptic strength on network activity amplitude and frequency are qualitatively equivalent in the CaSyn network which is indicated by symmetry of the heat plots (across the X=Y line) in Figure 4A,B. This symmetry results from the fact that the effective strength of synaptic interactions in the network is roughly proportional to a product of the synaptic strength and g-CAN. A transition from bursting to tonic spiking occurs when this effective excitation exceeds a certain critical value. This is why the bifurcation curve corresponding to a transition from rhythmic bursting to tonic spiking (a boundary between yellow and black in Figure 4A) looks like a hyperbola (g-CAN × synaptic strength=const).

Manipulations of synaptic strength (N·PSyn·12Wmax) and g-CAN have equivalent effects on network activity amplitude, frequency and recruitment of inspiratory neurons not involved in rhythm generation.

(A and B) Relationship between g-CAN  (mean values for the simulated populations), synaptic strength and the amplitude and frequency in the CaSyn network. Notice the symmetry about the X=Y line in panels A and B, which, indicates that changes in g-CAN and or synaptic strength are qualitatively equivalent. Synaptic strength was changed by varying Wmax. (C) Relationship between network activity amplitude and the reduction of g-CAN (blue) or synaptic strength (green). (D) Relationship between network frequency and the reduction of g-CAN (blue) or synaptic strength (green). (E and F) Decreasing g-CAN or synaptic strength de-recruits neurons by reducing the inspiratory drive potential, indicated by the amplitude of subthreshold depolarization, right traces. The solid blue and green lines in panels A and B represent the location in the 2D parameter space of the corresponding blue and green curves in C and D. The action potentials in the right traces of E and F are truncated to show the change in neuronal inspiratory drive potential. The parameters used for these simulations are CaSyng-Ca=[0,0], PCa=0.01, PSyn=1.0 and Wmax=var.

https://doi.org/10.7554/eLife.41555.005

We further investigated and compared the effect of reducing g-CAN or the synaptic strength on network activity amplitude and frequency as well as the effects on the recruitment of neurons not involved in rhythm generation (Figure 4C–F). To make this comparison, we picked a starting point in the 2D parameter space of g-CAN and synaptic strength where the network is bursting. Then in separate simulations, we gradually reduced either g-CAN or the synaptic strength to zero. We show that reducing either g-CAN or the synaptic strength have very similar effects on network activity amplitude and frequency (Figure 4C,D). Furthermore, de-recruitment of neurons in both cases is nearly identical (Figure 4E,F). Reducing either g-CAN or the synaptic strength decreases the excitatory input to the neurons during network oscillations which is a major component of the inspiratory drive potential. Therefore, in the CaSyn network, manipulations of g-CAN will affect the strength of the inspiratory drive potential in the rhythmic inspiratory neurons in a way that is equivalent to changing the synaptic strength of the network. In contrast, manipulations of g-CAN in the CaV network will only slightly affect the inspiratory drive potential due to changes in the average firing rate of active neurons (see Figure 3B).

Robustness of amplitude and frequency effects

We also examined if the effects are conserved in both the CaV and CaSyn networks over a range of network parameters. To test this, we investigated the dependence of network activity amplitude and frequency on g-CAN and average synaptic strength for CaSyn and CaV networks with high (PSyn=1) and low (PSyn=0.05) connection probabilities, and high (gCa=0.1 nS, PCa=0.1), medium (gCa=0.01 nS, PCa=0.01) and low (gCa=0.001 nS, PCa=0.005) strengths of calcium sources (Figure 5 and 6). We found that changing the synaptic connection probability and changing the strength of the calcium sources has no effect on the general relationship between g-CAN and the amplitude or frequency of bursts in the CaV or CaSyn networks. In other words, the general effect of increasing g-CAN on amplitude and frequency is conserved in both networks regardless of the synaptic connection probability or strength of the calcium sources. Increasing the strength of the calcium sources does, however, affects the range of possible g-CAN values where both networks produce rhythmic activity.

Robustness of amplitude and frequency effects to changes in g-CAN and synaptic strength in the CaV network for ‘high’ (left), ‘medium’ (middle) and ‘low’ (right) conductance of the voltage-gated calcium channel ICa as well as ‘high’ (top) and ‘low’ (bottom) network connection probabilities.

Amplitude and frequency are indicated by color (scale bar at right). Black regions indicate tonic network activity. Values of g-CAN indicated are the mean values for the simulated neuronal populations.

https://doi.org/10.7554/eLife.41555.006
Robustness of amplitude and frequency effects to changes in g-CAN and synaptic strength in the CaSyn network for ‘high’ (left), ‘medium’ (middle) and ‘low’ (right) calcium conductance in synaptic currents as well as ‘high’ (top) and ‘low’ (bottom) network connection probabilities.

Amplitude and frequency are indicated by color (scale bar at right). Black regions indicate tonic network activity. Values of g-CAN indicated are the mean values for the simulated neuronal populations.

https://doi.org/10.7554/eLife.41555.007

To summarize, in the CaV model, increasing g-CAN increases frequency, through increased excitability but has no effect on amplitude. In contrast, in the CaSyn model, increasing g-CAN slightly decreases frequency and increases amplitude. In this case, increasing g-CAN acts as a mechanism to increase the inspiratory drive potential and recruit previously silent neurons. Additionally, these features of the CaV and CaSyn models are robust and conserved across a wide range of network parameters.

Intracellular calcium transients activating ICAN primarily result from synaptically activated sources

In experiments where ICAN/TRPM4 was blocked by bath application of FFA or 9-phenanthrol (Koizumi et al., 2018) in vitro, the amplitude of network oscillations was strongly reduced and their frequency remained unchanged or was reduced insignificantly. Our model revealed that the effects of ICAN blockade on amplitude and frequency depend on the source(s) of intracellular calcium (see Figures 1 and 2). If the calcium influx is exclusively voltage-gated, our model predicts that ICAN blockade will have no effect on amplitude but reduce the frequency. In contrast, if the calcium source is exclusively synaptically gated, our model predicts that blocking ICAN will strongly reduce the amplitude and slightly increase the frequency. Therefore, a multi-fold decrease in amplitude, seen experimentally, is consistent with the synaptically driven calcium influx mechanisms, while nearly constant bursting frequency may be due to calcium influx through both voltage- and synaptically gated channels. Following the predictions above, to reproduce experimental data, we incorporated both mechanisms in the full model and inferred their individual contributions by finding the best fit. We found that the best match is observed (Figure 7) if synaptically mediated and voltage-gated calcium influxes comprise about 95% and 5% of the total calcium influx, respectively. We note that some experiments have shown a small (~20%) reduction of inspiratory burst frequency accompanying larger reductions in pre-BötC population activity amplitude with ICAN block (Peña et al., 2004). In our model, such perturbations of frequency can occur if the contribution of voltage-gated calcium influxes is larger than indicated above, or if neuronal background calcium concentration which partially activates ICAN is higher than specified in the model.

Experimental and simulated pharmacological blockade of ICAN by (A and C) 9-phenanthrol and (B and D) flufenamic acid (FFA).

Both voltage-gated and synaptic sources of intracellular calcium are included. Experimental blockade of ICAN (black) by 9-phenanthrol and FFA significantly reduce the (A and B) amplitude of network oscillations while having little effect on (C and D) frequency. The black line represents the mean and the gray band is the S.E.M. of experimental integrated XII output recorded from neonatal rat brainstem slices in vitro, reproduced from Koizumi et al., 2018. Simulated blockade of ICAN (red) closely matches the reduction in (A and B) amplitude of network oscillations and slight decrease in (C and D) frequency seen with 9-phenanthrol and FFA. Simulated and experimental blockade begins at the vertical dashed line. Blockade was simulated by exponential decay of g-CAN with the following parameters: 9-phenanthrol: ?Block=0.85tBlock=357s; FFA: ?Block=0.92tBlock=415s. The network parameters are: g-Ca=0.00175 (nS)PCa=0.0275, PSyn=0.05 and Wmax=0.096 (nS).

https://doi.org/10.7554/eLife.41555.008

INaP-dependent and Cai-ICAN-sensitive intrinsic bursting

In our model, we included INaP, ICAN as well as voltage-gated and synaptic mechanisms of Ca2+ influx. Activation of ICAN by CaSyn is the equivalent mechanism used in computational group-pacemaker models (Rubin et al., 2009aSong et al., 2015). Rhythmic burst generation and termination in our model, however, are dependent on INaP (Butera et al., 1999a). We investigated the sensitivity of intrinsic bursting in our model to INaP and calcium channel blockade (Figure 8). Intrinsic bursting was identified in neurons by zeroing the synaptic weights to simulate synaptic blockade. INaP and ICa blockade was simulated by setting g¯NaP and g¯Ca to 0nS. We found that after decoupling the network (Wmax=0) a subset of neurons remained rhythmically active (7% in this simulation) and that these were all neurons with a high INaP conductance. In these rhythmically active neurons, bursting was abolished in all neurons by INaP blockade. Interestingly, ICa blockade applied before INaP block abolished intrinsic bursting in three of the seven neurons and INaP block applied afterwards abolished intrinsic bursting in the remaining four neurons. Although only one rhythmogenic (INaP-based) mechanism exists in this model, bursting in a subset of these intrinsically bursting neurons is calcium sensitive, consistent with experimental observations of calcium-sensitive intrinsic bursters (Thoby-Brisson and Ramirez, 2001Del Negro et al., 2005; Peña et al., 2004). In calcium-sensitive bursters, Ca2+ blockade in our model abolishes bursting by reducing the intracellular calcium concentration and, hence, ICAN activation, which ultimately reduces excitability. We note that in our model the numbers of INaP-dependent and calcium-sensitive intrinsic bursters will vary depending on the mean and width of the INaP distribution and the background intracellular calcium concentration.

INaP-dependent and Ca2+-sensitive intrinsic bursting.

(A) From left to right, intrinsic bursters (pacemakers) are first identified by blocking synaptic connections. Cells whose activity is elminated under these conditions are non-pacemaker neurons. Then, calcium sensitive neurons are silenced and identified by ICa blockade. The remaining neurons are identifed as sensitive to INaP block. Top traces show the network output and Cell ID vs. g-NaP scatter plots that identify silent and bursting neurons under each condition. (B) INaP blockade after synaptic blockade eleminates bursting in all neurons. Therefore, all intrinsic bursters are INaP dependent. (C) Identification of calcium-sensitive and INaP-dependent as well as calcium-insensitive and INaP-dependent intrinsic bursters. Notice that only the neurons with the highest value of g-NaP are intrinsic bursters and that a subset of these neurons are sensitive to calcium blockade but all are dependent on INaP. The network parameters are: g-Ca=0.00175 (nS), PCa=0.0275, PSyn=0.05 and Wmax=0.096 (nS). The values of g-NaP given in the scatter plots indicate the magnitude of g-NaP for each neuron in the network and show the range of the g-NaP distribution.

https://doi.org/10.7554/eLife.41555.009

The rhythmogenic kernel

Our simulations have shown that the primary role of ICAN is amplitude but not oscillation frequency modulation with little or no effect on network activity frequency. Here we examined the neurons that remain active and maintain rhythm after ICAN blockade (Figure 9). We found that the neurons that remain active are primarily neurons with the highest g-NaP and that bursting in these neurons is dependent on INaP. Some variability exists and neurons with relatively low g-NaP value can remain active due to synaptic interactions while a neuron with a slightly higher g-NaP without sufficient synaptic input may become silent. These neurons that remain active after compete blockade of ICAN form a INaP-dependent kernel of a rhythm generating circuit.

ICAN blockade reveals an INaP-dependent rhythmogenic kernel.

The top traces show the network output at baseline, after ICAN blockade and INaP blockade. The bottom Cell ID vs. g-NaP scatter plots identify silent and bursting neurons in each conditon. Notice that only neurons with relitively high g-NaP remain active after ICAN block. The network parameters used are: g-Ca=0.00175 (nS), PCa=0.0275, PSyn=0.05 and Wmax=0.096 (nS). The values of g-NaP given in the scatter plots indicate the magnitude of g-NaP for each neuron in the network and show the range of the g-NaP distribution.

https://doi.org/10.7554/eLife.41555.010

Intracellular calcium transients and network activity during inhibition of ICAN/TRPM4

Dynamic calcium imaging has been used to assess activity of the population of pre-BöC excitatory neurons as well as individual pre-BötC neurons during pharmacological inhibition of ICAN/TRPM4 in vitro (Koizumi et al., 2018). These experiments indicate that network output activity amplitude and pre-BötC excitatory population-level intracellular calcium transients are highly correlated while the network oscillation frequency is not significantly perturbed. Interestingly, during ICAN/TRPM4 block, changes in calcium transients of individual neurons can differ significantly from the average population-level calcium transients. To assess if our model is consistent with these experimental results and gain additional insight into intracellular calcium dynamics during network activity, we analyzed simultaneous changes in the amplitude of network neuronal spiking activity, the average intracellular calcium concentration (Cai) of all network neurons, as well as Cai of individual neurons, with different network connection probabilities (PSyn) during simulated ICAN block (Figure 10). We found that regardless of PSyn, the network activity amplitude and average intracellular calcium concentration are highly correlated (Figure 10A,B). PSyn has no effect on the relationship between amplitude, calcium transients at the network level, or network oscillation frequency provided that the synaptic strength remains constant (NPSyn12Wmax=const). PSyn does, however, affect the change in the peak Cai in individual neurons. In a network with a high connection probability (PSyn=1) the synaptic current/calcium transient is nearly identical for all neurons and therefore the change in Cai during ICAN blockade is approximately the same for each neuron (Figure 10C). In a sparsely connected network, as proposed for the connectivity of the pre-BötC network (Carroll et al., 2013Carroll and Ramirez, 2013) the synaptic current and calcium influx are more variable and reflect the heterogeneity in spiking frequency of the pre-synaptic neurons (Figure 10). Interestingly, in a network with low connection probability (PSyn<0.1), the peak Cai transient in some neurons increases when ICAN is blocked (Figure 10E), consistent with our experimental results (Koizumi et al., 2018).

Changes of network activity amplitude, average network intracellular calcium concentration Cai amplitude, and single model neuron Cai amplitude during simulated ICAN blockade.

(A and B) Effect of ICAN block on network activity amplitude, network calcium amplitude and frequency for network connection probabilities (A) P = 1 and (B) P = 0.05. (C and D) Effect of ICAN block on changes in the magnitude of peak cellular calcium transients for network connection probabilities (C) PSyn=1 and (DPSyn=0.05. (E) Maximum, minimum and average change in the peak intracellular calcium transient of individual neurons as a function of synaptic connection probability. All curves in A through E are normalized to their baseline values. Synaptic weight was adjusted to keep the average synaptic strength (N·PSyn·12Wmax=const) constant. Notice that lowering the synaptic connection probability increases the variability in the peak intracellular calcium concentration during ICAN blockade. Interestingly, for connection probabilities below approximately 5%, blocking g-CAN can increase the peak calcium transient in a small subset of neurons. The network parameters used are: g-Ca=0.00175 (nS) and PCa=0.0275 and Wmax=var.

https://doi.org/10.7554/eLife.41555.011

Discussion

Establishing cellular and circuit mechanisms generating the rhythm and amplitude of respiratory oscillations in the mammalian brainstem pre-BötC has remained an unsolved problem of wide-spread interest in neurophysiology since this structure, essential for breathing to support mammalian life, was discovered nearly three decades ago (Smith et al., 1991). Our objective in this theoretical study was to re-examine and further define contributions of two of the main currently proposed neuronal biophysical mechanisms operating in pre-BötC excitatory circuits, specifically mechanisms involving ICAN activated by neuronal calcium fluxes (Thoby-Brisson and Ramirez, 2001; Peña et al., 2004; Del Negro et al., 2005; Mironov, 2008; Rubin et al., 2009a ) and voltage-dependent INaP in the circuit neurons (Butera et al., 1999aDel Negro et al., 2002; Koizumi and Smith, 2008; Koizumi and Smith, 2008). While these sodium- and calcium-based mechanisms have been studied extensively over the past two-decades and shown experimentally to be integrated in pre-BötC circuits, their actual roles in circuit operation are continuously debated and unresolved (Rybak et al., 2014). Both mechanisms have been proposed to be fundamentally involved in rhythm generation either separately or in combination, as plausibly shown from previous theoretical modeling studies (Butera et al., 1999aJasinski et al., 2013; Rubin et al., 2009a; Toporikova and Butera, 2011). Furthermore, the process of rhythm generation in pre-BötC circuits must be associated with an amplitude of excitatory circuit activity sufficient to drive downstream circuits to produce adequate respiratory motor output. Biophysical mechanisms involved in generating excitatory population activity amplitude have also not been established. Our analysis is motivated by the experimental observations obtained from neonatal rodent slices isolating pre-BötC circuits in vitro that inhibition of the endogenously active ICAN/TRPM4 strongly reduces the amplitude of network oscillations within pre-BötC circuits but has a small effect on oscillation frequency (Peña et al., 2004Koizumi et al., 2018). These findings challenge currently proposed ICAN-based models for rhythm generation in the isolated pre-BötC and indicate a functional organization of pre-BötC circuits, in terms of oscillatory frequency and amplitude generation, that needs to be defined.

We accordingly analyzed the role of ICAN and possible sources of intracellular calcium transients activating this conductance and found that the effect of simulated ICAN blockade on amplitude and frequency is highly dependent on the source(s) of intracellular calcium, which is also a central issue to be resolved. In the case where CaSyn is the primary intracellular calcium source, ICAN blockade generates a large reduction in network activity amplitude. In contrast, when CaV is the only intracellular calcium source, ICAN blockade has little effect on network activity amplitude and primarily affects the population bursting frequency that is caused by decreased excitability. Additionally, we show that activation of ICAN by CaSyn functions as a mechanism to augment the inspiratory drive potential and amplitude of population activity, and that this effect is similar to increasing the synaptic coupling strength within the network. Therefore, in the case of CaSyn, blockade of ICAN reduces the inspiratory drive potential causing de-recruitment of non-pacemaker rhythmic neurons and reduction of network activity amplitude. In a model where ICAN is activated by both CaV and CaSyn with contributions of 5% and 95% respectively, we show that simulated blockade of ICAN generates a large reduction in network population activity amplitude and a slight decrease in frequency. This closely reproduces experimental blockade of ICAN/TRPM4 by either 9-phenanthrol or FFA (Figure 7). Finally, we showed that the change in the peak calcium transients for individual neurons during ICAN blockade, particularly at relatively low network connection probabilities (PSyn∼<0.1), are consistent with experimental data.

Role of ICAN in the pre-BötC respiratory network

The hypothesis that ICAN is involved in generation of the inspiratory rhythm is based on experimental observations from in vitro mouse medullary slice preparations (Pace et al., 2007Mironov, 2008Del Negro et al., 2005; Del Negro et al., 2010; Peña et al., 2004; Thoby-Brisson and Ramirez, 2001), and in silico modeling studies (Jasinski et al., 2013Rybak et al., 2003; Toporikova and Butera, 2011). Theories of ICAN-dependent bursting rely on intracellular Ca2+ signaling mechanisms that have not been well defined.

Two models of ICAN-dependent rhythmic bursting in vitro have been proposed and are referred to as the 'dual pacemaker' and 'group pacemaker' models. In the dual pacemaker model, two types of pacemaker neurons are proposed that are either INaP-dependent (riluzole sensitive) or ICAN-dependent (Cd2+sensitive) intrinsic bursters (inspiratory 'pacemaker' neurons) (see Rybak et al., 2014 for review). In this model network, oscillations are thought to originate from these pacemaker neurons which through excitatory synaptic interactions synchronize bursting and drive activity of other rhythmic inspiratory neurons within the pre-BötC. Although pacemaker neurons sensitive to neuronal Ca2+ flux blockade through CaV have been reported (Peña et al., 2004; Thoby-Brisson and Ramirez, 2001, the source and mechanism driving intracellular Ca2+oscillations has not been definitely delineated. Computational models of ICAN-dependent pacemaker neurons rely on mechanisms for burst initiation and termination, for example IP3-dependent Ca2+oscillations (Toporikova and Butera, 2011; Del Negro et al., 2010), that have been questioned from recent negative experimental results (Beltran-Parrazal et al., 2012Toporikova et al., 2015). In favor of the dual pacemaker concept, experimental evidence has been presented that pharmacological block of both ICAN and INaP are necessary to disrupt rhythmogenesis under normoxic conditions in vitro (Peña et al., 2004). Other experimental results suggest that blocking INaP alone in vitro is sufficient (Koizumi and Smith, 2008; Toporikova et al., 2015), and we have discussed from a theoretical perspective how such discrepancies may depend on network excitability state and connectivity (Jasinski et al., 2013). Our present model defines more clearly parameters such as calcium flux sources and background calcium levels that can influence the relative contributions of ICAN in various network states to be further explored in model simulation studies.

In versions of the 'group pacemaker' model (Rubin et al., 2009aDel Negro et al., 2010; Feldman and Del Negro, 2006), network oscillations are initiated through recurrent synaptic excitation that trigger postsynaptic Ca2+ influx. Subsequent ICAN activation generates membrane depolarization (inspiratory drive potential) to drive neuronal bursting. Synaptically triggered Ca2+ influx and the contribution of ICAN to the inspiratory drive potential of individual pre-BötC neurons are experimentally supported (Mironov, 2008Pace et al., 2007); however, the mechanism of burst termination remains unclear. Again, the computational group-pacemaker models that have been explored (Rubin et al., 2009a) rely on speculative mechanisms for burst termination that need to be tested experimentally, and in some cases lack key biophysical features of the pre-BötC neurons such as voltage-dependent frequency control and expression of INaP.

In our model, we showed that blockade of either ICAN or synaptic interactions produce qualitatively equivalent effects on network population activity amplitude and frequency when the calcium transients are primarily generated from synaptic sources (Figure 4). Consequently, our model predicts that blockade of ICAN or synaptic interactions in the isolated pre-BötC in vitro will produce comparable effects on amplitude and frequency. This is the case as Johnson et al. (1994) showed that gradual blockade of synaptic interactions by low calcium solution significantly decreases network activity amplitude while having little effect on frequency, similar to the experiments where the ICAN channel TRPM4 is blocked with 9-phenanthrol (Koizumi et al., 2018). We note that complete blockade of ICAN in our model can ultimately abolish synchronized network oscillations due to weakened excitatory synaptic transmission, which results in neuronal de-recruitment and desynchronization of the network, particularly when synaptic strength is low (Figure 4). Thus, ICAN plays a critical role in network activity synchronization that determines the ability of the pre-BötC excitatory network to produce rhythmic output, depending on synaptic strength.

Overall, our new model simulations for the isolated pre-BötC excitatory network suggest that the role of ICAN/TRPM4 activation is to amplify excitatory synaptic drive in generating the amplitude of inspiratory population activity, essentially independent of the biophysical mechanism generating inspiratory rhythm. We note that the recent experiments have also shown that in the more intact brainstem respiratory network that ordinarily generates patterns of inspiratory and expiratory activity, endogenous activation of ICAN/TRPM4 appears to augment the amplitude of both inspiratory and expiratory population activity, and hence these channels are fundamentally involved in inspiratory-expiratory pattern formation (Koizumi et al., 2018).

Intracellular calcium dynamics and network activity during inhibition of ICAN/TRPM4

We analyzed the correlation between calcium transients and inspiratory activity of individual inspiratory neurons as well as the entire network, particularly since dynamic calcium imaging has been utilized to assess activity of individual cells and populations of pre-BötC excitatory neurons in vitro during pharmacological inhibition of ICAN/TRPM4 (Koizumi et al., 2018). We show that intracellular calcium transients are highly correlated with network and cellular activity across the duration of an ICAN blockade simulation, consistent with experimental observations.

Additionally, we examined the relative change in the peak calcium transients in single neurons as a function of ICAN conductance. We show that in a subset of neurons the peak calcium transient increases with reduced ICAN. This result is surprising but is also supported by the calcium imaging data (Koizumi et al., 2018). This occurs in neurons that receive most of their synaptic input from pacemaker neurons and our analyses suggest this is possible in sparse networks, that is with low connection probability. In pacemaker neurons, ICAN blockade leads to a reduction of their excitability resulting in an increased value of INaP inactivation gating variable at the burst onset. Thus, during the burst, the peak action potential frequency and the synaptic output from these neurons is increased with ICAN blockade. Consequently, neurons that receive synaptic input from pacemaker neurons will see an increase in their peak calcium transients. In most neurons, however, synaptic input is received primarily from non-pacemaker rhythmic neurons. Since ICAN blockade de-recruits these non-pacemaker neurons, the synaptic input and subsequent calcium influx in most of these cells decreases. Therefore, our model predicts that in a sparse network, which has been proposed for pre-BötC network connectivity (Carroll et al., 2013Carroll and Ramirez, 2013), blocking ICAN results in very diverse responses at the cellular level with an overall tendency to reduce intracellular calcium transients such that the amplitude of these transients averaged over the entire population decreases during ICAN blockade, while their burst frequency is essentially unchanged, as found from the calcium imaging experiments (Koizumi et al., 2018).

Synaptic calcium sources

Our model suggests that calcium transients in the pre-BötC are coupled to excitatory synaptic input, that is pre-synaptic glutamate release and binding to post-synaptic glutamate receptors triggers calcium entry. The specific mechanisms behind this process are unclear; however, this is likely dependent on specific types of ionotropic or metabotropic glutamate receptors.

There are three subtypes of ionotropic glutamate receptors, N-methyl-D-aspartate (NMDA), Kainate (KAR), and a-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA), all of which are expressed in the pre-BötC (Paarmann et al., 2000) and have varying degrees of calcium permeability. NMDA and AMPA are unlikely candidates for direct involvement in synaptically mediated calcium influx in the pre-BötC. Pharmacological blockade of NMDA receptors does not consistently effect the amplitude or frequency of XII motor output (Lieske and Ramirez, 2006a; Morgado-Valle and Feldman, 2007; Pace et al., 2007) and AMPA receptors in the pre-BötC show high expression of the subunit GluR2, which renders the AMPA ion channel pore impermeable to Ca2+ (Paarmann et al., 2000). It is possible, however, that AMPA mediated depolarization may trigger calcium influx indirectly through the voltage-gated calcium channel activation on the post-synaptic terminal. The latter may contribute to synaptically triggered calcium influx as blockade of P/Q-type (but not L- N-type) calcium channels reduces XII motor output driven from the pre-BötC in the in vitro mouse slice preparations from normal animals (Lieske and Ramirez, 2006a; Koch et al., 2013).

Calcium permeability through KAR receptors is dependent on subunit expression. The KAR subunit GluK3 is highly expressed in the pre-BötC (Paarmann et al., 2000) and is calcium permeable (Perrais et al., 2009) making it a possible candidate for synaptically mediated calcium entry. Furthermore, GluK3 is insensitive to tonic glutamate release and only activated by large glutamate transients (Perrais et al., 2009). Consequently, GluK3 may only be activated when receiving synaptic input from a bursting presynaptic neuron which would presumably generate large glutamate transients. The role of GluK3 in the pre-BötC has not been investigated.

Metabotropic glutamate receptors (mGluR) indirectly activate ion channels through G-protein mediated signaling cascades. Group 1 mGluRs which include mGluR1 and mGluR5 are typically located on post-synaptic terminals (Shigemoto et al., 1997) and activation of group 1 mGluRs is commonly associated with calcium influx through calcium permeable channels (Berg et al., 2007; Endoh, 2004; Mironov, 2008) and calcium release from intracellular calcium stores (Pace et al., 2007).

In the pre-BötC, mGluR1/5 are thought to contribute to calcium influx by triggering the release of calcium from intracellular stores (Pace et al., 2007) and/or the activation of the transient receptor potential C3 (TRPC3) channel (Ben-Mabrouk and Tryba, 2010). Blockade of mGluR1/5 reduces the inspiratory drive potential in pre-BötC neurons (Pace et al., 2007) without significant perturbation of inspiratory frequency (Pace et al., 2007Lieske and Ramirez, 2006a), which is consistent with the effects of ICAN/TRPM4 blockade (Koizumi et al., 2018). TRPC3 is a calcium permeable channel (Thebault et al., 2005) that is associated with calcium signaling (Hartmann et al., 2011), store-operated calcium entry (Kwan et al., 2004), and synaptic transmission (Hartmann et al., 2011). TRPC3 is activated by diacylglycerol (DAG) (Clapham, 2003), which is formed after synaptic activation of mGluR1/5. TRPC3, which is highly expressed in pre-BötC glutamatergic inspiratory neurons, is often co-expressed with TRPM4 (Koizumi et al., 2018) and was hypothesized to underlie ICAN activation in the pre-BötC (Ben-Mabrouk and Tryba, 2010) and other brain regions (Amaral and Pozzo-Miller, 2007Zitt et al., 1997). Furthermore, TRPC3 and ICAN have been shown to underlie slow excitatory post-synaptic current (sEPSC) (Hartmann et al., 2008Hartmann et al., 2011). This is consistent with our model since ICAN activation is dependent on synaptically triggered calcium entry, and the calcium dynamics are slower than the fast AMPA based current ISyn. Therefore, in our model, ICAN decays relatively slowly and, hence, can be treated as a sEPSC.

The effect of TRPC3 blockade by 3-pyrazole on pre-BötC network activity amplitude is remarkably similar to that with blockade of TRPM4 (Koizumi et al., 2018). This suggests that the ICAN/TRPM4 activation may be dependent on/coupled to TRPC3. A possible explanation is that TRPC3, which is calcium permeable, mediates synaptically triggered calcium entry. It is also likely that TRPC3 plays a role in maintaining background calcium concentration levels. We tested this hypothesis by simulating the blockade of synaptically-triggered calcium influx while simultaneously lowering the background calcium concentration (Figure 11). These simulations generated large reductions in activity amplitude with no effect on frequency which are consistent with data from experiments where TRPC3 is blocked using 3-pyrazole (Koizumi et al., 2018). This indirectly suggests that TRPC3 may be critical for synaptically-triggered calcium entry and subsequent ICAN activation.

Comparison of experimental (black) and simulated (red) TRPC3 blockade (by CaSyn block) on network activity amplitude (A) and frequency (B).

Simulated and experimental blockade begins at the vertical dashed line. The black line represents the mean and the gray band represents the S.E.M. of experimental integrated XII output recorded from neonatal rat brainstem slices in vitro, reproduced from Koizumi et al., 2018. Blockade was simulated by exponential decay of PCa with the following parameters: 3-pyrazole: ?Block=1.0τBlock=522.5s. The network parameters are: g-Ca=0.00175 (nS), PCa=0.0275, PSyn=0.05 and Wmax=0.096 (nS).

https://doi.org/10.7554/eLife.41555.012

INaP -dependent rhythmogenic kernel

INaP is a conductance present ubiquitously in pre-BötC inspiratory neurons, and is established to underlie intrinsic oscillatory neuronal bursting in the absence of excitatory synaptic interactions in neurons with sufficiently high INaP conductance densities (Del Negro et al., 2002Koizumi and Smith, 2008; Koizumi and Smith, 2008; Yamanishi et al., 2018). Accordingly, we randomly incorporated this conductance in our model excitatory neurons from a uniform statistical distribution to produce heterogeneity in INaP conductance density across the population. Our simulations indicate that the circuit neurons mostly with relatively high INaP conductance values underlie rhythm generation and remain active after compete blockade of ICAN in our model network, thus forming a INaP-dependent rhythmogenic kernel, including some neurons with intrinsic oscillatory bursting behavior when synaptically uncoupled.

As noted above, the rhythmogenic properties of individual neurons depend on whether their INaP conductance is high enough and, therefore, the number of intrinsic bursters in the model is defined by the width of this conductance distribution over the population. However, the critical value of INaP conductance and intrinsic bursting properties in general were also shown theoretically and experimentally to critically depend on the conductances of other ionic channels (e.g. leak and delayed rectifier potassium conductances) and extracellular ion concentrations such as potassium concentration (Bacak et al., 2016bKoizumi and Smith, 2008; Rybak et al., 2003). Therefore, we believe that in reality the specific composition of the rhythmogenic kernel and its oscillatory capabilities strongly depend on the existing combination of neuronal conductances and on the in vitro experimental conditions.

Recently, it has become apparent that there is functional heterogeneity within pre-BötC excitatory circuits, including distinct subpopulations of neurons involved in generating periodic sighs (Toporikova et al., 2015; Li et al., 2016), arousal (Yackle et al., 2017), and the subpopulations generating regular inspiratory activity. Activity of the normal inspiratory and sigh-generating subpopulations in the pre-BötC isolated in vitro is proposed to be dependent on activation of INaP (Toporikova et al., 2015). Our experimental and modeling results suggest that within the normal inspiratory population, there are subpopulations distinguished by their role in rhythm versus amplitude generation due to biophysical properties: there is a ICAN/TRPM4-dependent recruitable population of excitatory neurons for burst amplitude generation and the INaP-dependent rhythmogenic kernel population. The spatial arrangements of these two synaptically interconnected excitatory populations within the pre-BötC are currently unknown, and it remains an important experimental problem to identify the cells constituting the rhythmogenic kernel and their biophysical properties. This should now be possible, since our analysis and experimental results suggest that the rhythmically active neurons of the kernel population can be revealed and studied after pharmacologically inhibiting the ICAN/TRPM4-dependent inspiratory burst-generating population.

A 'burstlet theory' for emergent network rhythms has recently been proposed to account for inspiratory rhythm and pattern generation in the isolated pre-BötC in vitro (Kam et al., 2013Del Negro et al., 2018). This theory postulates that a subpopulation of excitatory neurons generating small amplitude oscillations (burstlets) functions as the inspiratory rhythm generator that drives neurons that generate the larger amplitude, synchronized inspiratory population bursts. This concept emphasizes that subthreshold neuronal membrane oscillations need to be considered and that there is a neuronal subpopulation that functions to independently form the main inspiratory bursts. This is similar to our concept of distinct excitatory subpopulations generating the rhythm versus the amplitude of inspiratory oscillations. Biophysical mechanisms generating rhythmic burstlets and the large amplitude inspiratory population bursts in the burstlet theory are unknown, and the general problem of understanding the dynamic interplay of circuit interactions and cellular biophysical processes in the generation of population-level bursting activity has been highlighted (Richter and Smith, 2014Ramirez and Baertsch, 2018b). We have identified a major Ca2+-dependent conductance mechanism for inspiratory burst amplitude (pattern) generation and show theoretically how this mechanism may be coupled to excitatory synaptic interactions and is independent of the rhythm-generating mechanism. We also note that a basic property of INaP is its ability to generate subthreshold oscillations and promote burst synchronization (Butera et al., 1999b; Bacak et al., 2016a). However, in contrast to our proposal for the mechanisms operating in the kernel rhythm-generating subpopulation, INaP with its favorable voltage-dependent and kinetic autorhythmic properties– is not proposed to be a basic biophysical mechanism for rhythm generation in the burstlet theory (Del Negro et al., 2018).

We emphasize that the above discussions regarding the role of INaP pertain to the excitatory circuits in the isolated pre-BötC including in more mature rodent experimental preparations in situ where inspiratory rhythm generation has also been shown to be dependent on INaP (Smith et al., 2007). The analysis is more complex when the pre-BötC is embedded within interacting respiratory circuits in the intact nervous system generating the full complement of inspiratory and expiratory phase activity (Lindsey et al., 2012Ramirez and Baertsch, 2018a; Richter and Smith, 2014), where rhythmogenesis is tightly controlled by inhibitory circuit interactions, including via local inhibitory circuits in the pre-BötC (Harris et al., 2017Ausborn et al., 2018; Baertsch et al., 2018; Ramirez and Baertsch, 2018a), and the contribution of INaP kinetic properties alone in setting the timing of inspiratory oscillations is diminished (Smith et al., 2007Richter and Smith, 2014; Rubin et al., 2009b). Extending our analysis to consider inhibitory circuit interactions with the excitatory subpopulations generating oscillation frequency and amplitude that we propose will provide additional insight into biophysical mechanisms controlling these two processes.

Conclusions

Based on our computational model, distinct biophysical mechanisms are involved in generating the rhythm and amplitude of inspiratory oscillations in the isolated pre-BötC excitatory circuits. According to this model, inspiratory rhythm generation arises from a group of INaP-dependent excitatory neurons, including cells with intrinsic oscillatory bursting properties, that form a rhythmogenic kernel. Rhythmic synaptic drive from these neurons triggers post-synaptic calcium transients, ICAN activation, and subsequent membrane depolarization which drives rhythmic bursting in the rest of the population of inspiratory neurons. We showed that activation of ICAN by synaptically-driven calcium influx functions as a mechanism that amplifies the excitatory synaptic input to generate the inspiratory drive potential and population activity amplitude in these non-rhythmogenic neurons. Consequently, reduction of ICAN causes a robust decrease in overall network activity amplitude via de-recruitment of these burst amplitude-generating neurons without substantial perturbations of the inspiratory rhythm. Thus, ICAN plays a critical role in generating the amplitude of rhythmic population activity, which is consistent with the results from experimental inhibition of ICAN/TRPM4 channels (Peña et al., 2004Koizumi et al., 2018). Our model provides a theoretical explanation for these experimental results and new insights into the biophysical operation of pre-BötC excitatory circuits. The theoretical framework that we have developed here should provide the bases for further exploration of biophysical mechanisms operating in the mammalian respiratory oscillator.

Materials and methods

Model description

The model describes a network of N=100 synaptically coupled excitatory neurons. Simulated neurons are comprised of a single compartment described using a Hodgkin Huxley formalism. For each neuron, the membrane potential Vm is given by the following current balance equation:

CmdVmdt+INa+IK+ILeak+INaP+ICAN+ICa+ISyn=0

where Cm is the membrane capacitance, INa, IK, ILeak, INaP, ICANICa and ISyn are ionic currents through sodium, potassium, leak, persistent sodium, calcium activated non-selective cation, voltage-gated calcium, and synaptic channels, respectively. Description of these currents, synaptic interactions, and parameter values are taken from Jasinski et al. (2013). The channel currents are defined as follows:

INa=g¯Na·mNa3·hNa·(Vm-ENa)
IK=g¯K·mK4·(Vm-EK)
ILeak=g¯Leak(VmELeak)
INaP=g¯NaP·mNaP·hNaP·(Vm-ENa)
ICAN=g¯CAN·mCAN·(Vm-ECAN)
ICa=g¯Ca·mCa·hCa·(Vm-ECa)
ISyn=gSyn·Vm-ESyn

where g¯i is the maximum conductance, Ei is the reversal potential, mi and hi are voltage dependent gating variables for channel activation and inactivation, respectively, and i?{Na,K,Leak,NaP,CAN,Ca,Syn}. The parameters g¯i and Ei are given in Table 1.

For INa, IK, INaP, and ICa, the dynamics of voltage-dependent gating variables mi, and hi are defined by the following differential equation:

t?V·d?dt=?8V-?; ??{mi,hi}

where steady state activation/inactivation ?8 and time constant t? are given by:

?8V=1+e-(V-V?1/2)/k?-1
t?(V)=t?max/cosh ?((V-Vt?1/2)/kt?).

For the voltage-gated potassium channel, steady state activation mK8(V) and time constant tmK(V) are given by:

mK8V=a8Va8V+ß8V
tmK(V)=1/(a8(V)+ß8(V))

where

a8(V)=Aa·(V+Ba)/(1-exp(-(V+Ba)/?a))
ß8(V)=Aß·exp(-(V+Bß)/?ß).

The parameters V?1/2, Vt?1/2, ??, κτητηmax, Aa, Aß, Ba, Bß, ?a, and ?ß are given in Table 1. ICAN activation is dependent on the intracellular calcium concentration Cain and is given by:

mCAN=1/(1+(Ca1/2/Cain)n).
Table 1
Model parameter values.

The channel kinetics, intracellular Ca2+ dynamics and the corresponding parameter values, were derived from previous models (see Jasinski et al., 2013) and the references therein).

https://doi.org/10.7554/eLife.41555.013
ChannelParameters
INag-Na=150.0 nS, ENa=55.0 mV,Vm1/2=43.8mV, km=6.0 mV,Vτm1/2=-43.8 mV, kτm=14.0 mV, τmmax=0.25 ms,Vh1/2=67.5mV, kh=-10.8 mV,Vτh1/2=-67.5 mV, kτh=12.8 mV, τhmax=8.46 ms
IKg-K=160.0 nS, EK=-94.0 mV,
Aα=0.01, Bα=44.0 mV, κα=5.0 mV
Aβ=0.17, Bβ=49.0 mV, κβ=40.0 mV
ILeak

g-Leak=2.5 nS, ELeak=-68.0 mV

INaPgNaP[0.0,5.0]nS,
Vm1/2=47.1mV, km=3.1 mV,
Vτm1/2=-47.1 mV, kτm=6.2 mV, τmmax=1.0 ms,
Vh1/2=60.0mV, kh=-9.0 mV,
Vτh1/2=-60.0 mV, kτh=9.0 mV, τhmax=5000 ms
ICANgCAN[0.5,1.5]nS, ECAN=0.0 mV,
Ca1/2=0.00074 mM, n=0.97
ICag-Ca=0.01 nS, ECa=RT/FlnCaout/Cain,
R=8.314 J/(molK), T=308.0 K,
F=96.485 kC/mol, Caout=4.0 mM
Vm1/2=27.5mV, km=5.7 mV, τm=0.5 ms,
Vh1/2=52.4mV, kh=-5.2 mV, τh=18.0 ms
CainαCa=2.510-5mM/fCPCa=0.01, Camin= 1.010-10mM, τCa=50.0ms
ISyngTonic=0.31 nS, ESyn=-10.0 mV, τSyn=5.0 ms

The parameters Ca1/2 and n, given in Table 1, represent the half-activation calcium concentration and the Hill Coefficient, respectively.

Calcium enters the neurons through voltage-gated calcium channels (CaV) and/or synaptic channels (CaSyn), where a percentage (PCa) of the synaptic current (ISyn) is assumed to consist of Ca2+ ions. A calcium pump removes excess calcium with a time constant tCa and sets the minimum calcium concentration Camin. The dynamics of Cain is given by the following differential equation:

dCaindt=-aCa(ICa+PCa·Isyn)-(Cain-Camin)/tCa.

The parameters aCa is a conversion factor relating current and rate of change in Cain, see Table 1 for parameter values.

The synaptic conductance of the ith neuron (gSyni) in the population is described by the following equation:

gSyni=gTonic+?j,nwji·Cji·H(t-tj,n)·e-(t-tj,n)/tsyn

where wji is the weight of the synaptic connection from cell j to cell i, C is a connectivity matrix (Cji=1 if neuron j makes a synapse on neuron i, and Cji=0 otherwise), H(.) is the Heaviside step function, t is time, tSyn is the exponential decay constant and tj,n is the time at which an action potential n is generated in neuron j and reaches neuron i.

To account for heterogeneity of neuron properties within the network, the persistent sodium current conductance, g-NaP, for each neuron was assigned randomly based on a uniform distribution over the range 0.0,5.0 nS which is consistent with experimental measurements (Rybak et al., 2003; Koizumi and Smith, 2008; Koizumi and Smith, 2008). We also uniformly distributed g-CAN over the range [0.5,1.5] nS , however, simulation results did not depend on whether we used such a distribution, or assigned g-CAN for all neurons to the same value of 1.0 nS, which is the mean of this distribution. In simulations where g-CAN was varied, we multiplied g-CAN for each neuron by the same factor. This factor was used as a control parameter for all such simulations, and shown as a percentage of the baseline g-CAN or as the mean gCAN values for the population in figures. The weight of each synaptic connection was uniformly distributed over the range wji?[0,Wmax] where Wmax ranged from 0.0 to 1.0 nS depending on the network connectivity and specific simulation. The elements of the network connectivity matrix, Cji, are randomly assigned values of 0 or 1 such that the probability of any connection between neuron j and neuron i being 1 is equal to the network connection probability PSyn. We varied the connection probability over the range PSyn[0.05,1.0], however, a value of PSyn=0.05 was used in most simulations.

Data analysis and definitions

The time of an action potential was defined as when the membrane potential of a neuron crosses -35mV in a positive direction. The network activity amplitude and frequency were determined by identifying peaks and calculating the inverse of the interpeak interval in histograms of network spiking. Network histograms of the population activity were calculated as the number of action potentials generated by all neurons per 50 ms bin per neuron with units of spikes/s. The number of recruited neurons is defined as the peak number of neurons that spiked at least once per bin during a network burst. The average spike frequency of recruited neurons is defined as the number of action potentials per bin per recruited neuron with units of spikes/s. The average network resting membrane potential was defined as the average minimum value of Vm in a 500 ms window following a network burst. The average inactivation of the persistent sodium current at the start of each burst was defined by the maximum of the average value of hNaP in a 500-ms window before the peak of each network burst. The average inactivation of the persistent sodium current at the end of each burst was defined by the maximum of the average value of hNaP in a 500-ms window after the peak of each network burst. Synaptic strength is defined as the number of neurons in the network multiplied by the connection probability multiplied by the average weight of synaptic connections (N·PSyn·12Wmax). Pacemaker neurons were defined as neurons that continue bursting intrinsically after complete synaptic blockade. Follower neurons were defined as neurons that become silent after complete synaptic blockade. The inspiratory drive potential is defined as the envelope of depolarization that occurs in neurons during the inspiratory phase of the network oscillations (Morgado-Valle et al., 2008).

Characterization of ICAN in regulating network activity amplitude and frequency in CaV and CaSyn models

To characterize the role of ICAN in regulation of network activity amplitude and frequency we slowly increased the conductance (g-CAN) in our simulations from zero until the network transitioned from a rhythmic bursting to a tonic (non-bursting) firing regime. To ensure that the effect(s) are robust, these simulations were repeated over a wide range of synaptic weights, synaptic connection probabilities, and strengths of the intracellular calcium transients from CaV or CaSyn sources. Changes in network activity amplitude were further examined by plotting the number of recruited neurons and the average action potential frequency of recruited neurons versus g-CAN.

Simulated pharmacological manipulations

In simulations that are compared with experimental data, both CaV and CaSyn calcium sources are included. Pharmacological blockade of ICAN was simulated by varying the conductance, g-CAN according to a decaying exponential function

gCAN(t)=gCANmaxγblock(1et/τblock).

The percent block ?block, decay constant tblock and the maximum ICAN conductance gCANmax were adjusted to match the experimental changes in network amplitude. The synaptic weight of the network was chosen such that at g-CAN=0 the network activity amplitude was close to 20% of maximum. To reduce the computational time, the duration of ICAN block simulations was one tenth of the total of experimental durations. For comparison, the plots of normalized change in amplitude and frequency of the simulations were stretched over the same time-period as experimental data. Increasing the simulation time had no effect on our results (data not shown).

Comparison with calcium imaging data

To allow comparisons with network and cellular calcium imaging data, we analyzed rhythmic calcium transients from our simulations. Single cell calcium signals are represented by [Ca]i. The network calcium signal was calculated as the average intracellular calcium concentration in the network (1N[Ca]i/N).

Integration methods

All simulations were performed locally on an 8-core Linux-based operating system or on the high-performance computing cluster Biowulf at the National Institutes of Health. Simulation software was custom written in C++. Numerical integration was performed using the exponential Euler method with a fixed step-size (Δt) of 0.025ms. In all simulations, the first 50 s of simulation time was discarded to allow for the decay of any initial condition-dependent transients.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
    Mechanisms leading to rhythm cessation in the respiratory PreBötzinger complex due to piecewise cumulative neuronal deletions
    1. H Song
    2. JA Hayes
    3. NC Vann
    4. M Drew LaMar
    5. CA Del Negro
    (2015)
    eNeuro, 2, 10.1523/ENEURO.0031-15.2015, 26465010.
  63. 63
  64. 64
  65. 65
    Sigh and eupnea rhythmogenesis involve distinct interconnected subpopulations: a combined computational and experimental study
    1. N Toporikova
    2. M Chevalier
    3. M Thoby-Brisson
    (2015)
    eNeuro, 2, 10.1523/ENEURO.0074-14.2015, 26464980.
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70

Decision letter

  1. Ronald L Calabrese
    Senior Editor; Emory University, United States
  2. Jan-Marino Ramirez
    Reviewing Editor; Seattle Children's Research Institute and University of Washington, United States
  3. Jan-Marino Ramirez
    Reviewer; Seattle Children's Research Institute and University of Washington, United States
  4. Guillaume Drion
    Reviewer

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Biophysical mechanisms in the mammalian respiratory oscillator re-examined with a new data-driven computational model" for consideration by eLife. Your article has been reviewed by two peer reviewers, including Jan-Marino Ramirez as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Ronald Calabrese as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Guillaume Drion (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

The manuscript by Phillips et al. explores the interplay between two inward currents (ICAN and INaP) that have been implicated in respiratory rhythm generation within the pre-Bötzinger complex. Using computational modeling the authors carefully characterize the relative rhythmogenic contributions of each of these currents, which are dynamic, can differ in individual neurons and can change when one current is blocked, or when the intracellular calcium concentration changes. The authors use a computational model made of 100 synaptically coupled excitatory neurons with random connectivity and random INaP conductance (and ICAN conductance?). Each neuron model is conductance-based and includes 5 voltage-gated conductances, 1 leak conductance, intracellular calcium dynamics and synaptic currents. The paper compares two different model configurations: a configuration where voltage-gated calcium channels are the main source of calcium for ICAN activation (CaV network) and a configuration where ICAN is activated by calcium coming from synaptically activated sources (CaSyn network). First, analyses on the role of ICAN conductance on network amplitude and frequency in both configurations and comparison with recent experimental data suggest that ICAN activation mostly relies on synaptically-activated calcium sources, which confirms previous studies. Second, the paper shows that INaP is critical for network rhythmic activity: it imbues some neurons of the network with spontaneous bursting capabilities, forming an INaP-dependent rhythmogenic kernel. We find this study important, in part because the dynamic nature could also explain why there are so many discrepancies in the field regarding the role of these two inward currents.

Title:

In the title, the paper mentions the use of a "data-driven computational model". I was therefore expecting to see a model whose parameters were set according to experimental data, which is not the case in this paper (or at least not a contribution of this paper). I understand that the model structure was designed to investigate the mechanisms that could explain recent experimental data, which I appreciate, but calling the model "data-driven" seems misleading to me (even if data are used to quantify the contributions of CaV and CaSyn in Figure 7).

Essential revisions:

1) The authors focused only on the fact that blocking ICAN has a small frequency effect but large amplitude effect, which is something that has been previously reported in 2004 (Pena et al., 2004), and which these researchers reproduced in 2018. The authors don't address the finding that blocking INaP alone has also no effect on the frequency (Pena et al., 2004), while others show that blocking INaP alone can block rhythmogenesis under control conditions. By contrast Pena et al. found that this is only the case in severe hypoxia. The modeling could explain these different experimental findings. For example the intracellular calcium concentrations or the network state could differ between different labs. Unfortunately, the authors don't pitch the paper like this, which seems like a missed opportunity. Thus, the authors should try to be less polarizing and report data as is (and as detailed below in more detail). There is no need to try to push the agenda that only INaP is rhythmogenic when we know that ICAN has a big role in rhythmogenesis by regulating synaptic excitatory transmission.

2) The authors should add a discussion regarding the neuron model and the robustness of the network activity. First, the model generates bursts of action potentials having an amplitude of around 30-40mV, as shown in Figure 4 and Figure 10. But in the experimental data of Koizumi et al., the amplitude of action potentials is more about roughly 80-100 mV (Figure 3 of the paper). Why is there such a big difference between the two? The spike height is indeed often used as a major criterion to assess the validity of a computational model at the cellular level. We agree that it is not a critical feature of the model used in this paper, but the issue deserves to be discussed. In particular, what happens to the bursting pattern if you attempt to reach a physiologically plausible spike height by increasing the transient sodium conductance (which might also require an increase in delayed-rectifier potassium conductance)? We suspect that bursting would be very fragile to an increase in Na conductance, due to the specific dynamical mechanisms at play in the model.

3) Along the same lines, Figure 4A suggests that the region of network bursting is quite narrow as compared to the region of network tonic spiking (especially if one increases the parameter range). What is the robustness of this bursting region against variability/heterogeneity in the other conductance values? Do you think that this result implies that physiological mechanisms need to tightly regulate the values of the conductances to maintain a viable network bursting activity (for instance a specific balance between synaptic strength and ICAN conductance), or could it be possible that this fragility is an artifact of the computational model?

4) The authors should succinctly describe the model at the beginning of the Results section. Indeed, it was difficult for us to understand the configuration of the model without having to go through the Materials and methods section, which is at the end of the manuscript and very detailed. It is important that the reader has a clear view of how neurons are modeled (ion channels), how many neurons constitute the circuit, how you compute "circuit activity amplitude", etc. before going through the analyses. Please also make sure to define every specific term when it is introduced. For instance, when describing how network probability is computed, Wmax is not defined, and it cannot be easily related to the "average weight of synaptic connections" without again going through the Materials and methods section. Also, you mention in the Abstract that ICAN conductance is distributed randomly, but not in the Materials and methods section. Which one is correct?

https://doi.org/10.7554/eLife.41555.017

Author response

Title:

In the title, the paper mentions the use of a "data-driven computational model". I was therefore expecting to see a model whose parameters were set according to experimental data, which is not the case in this paper (or at least not a contribution of this paper). I understand that the model structure was designed to investigate the mechanisms that could explain recent experimental data, which I appreciate, but calling the model "data-driven" seems misleading to me (even if data are used to quantify the contributions of CaV and CaSyn in Figure 7).

We do not see the problem for the reasons that the reviewers state and do not think that calling the model “data-driven” is misleading. The modeling was entirely motivated by new data requiring a theoretical and conceptual explanation including model fits to some of the key data.

Essential revisions:

1) The authors focused only on the fact that blocking ICAN has a small frequency effect but large amplitude effect, which is something that has been previously reported in 2004 (Pena et al., 2004), and which these researchers reproduced in 2018. The authors don't address the finding that blocking INaP alone has also no effect on the frequency (Pena et al., 2004), while others show that blocking INaP alone can block rhythmogenesis under control conditions. By contrast Pena et al. found that this is only the case in severe hypoxia. The modeling could explain these different experimental findings. For example the intracellular calcium concentrations or the network state could differ between different labs. Unfortunately, the authors don't pitch the paper like this, which seems like a missed opportunity. Thus, the authors should try to be less polarizing and report data as is (and as detailed below in more detail). There is no need to try to push the agenda that only INaP is rhythmogenic when we know that ICAN has a big role in rhythmogenesis by regulating synaptic excitatory transmission.

We focused on independence of the amplitude and the frequency of oscillations which is mathematically counterintuitive and represents a major qualitative effect observed experimentally that must be indicating something important about the dynamical/biophysical organization and operation of the isolated pre-BötC network. This problem has not been previously addressed. Reference to Pena et al. (2004) was added in the Introduction and throughout. While amplitude modulation associated with experimental attenuation of ICAN conductance has been reported, these and other studies referenced have emphasized the potential role of ICAN in intrinsic neuronal bursting and rhythm generation by pre-BötC circuits, justifiably since calcium-dependent/ICAN intrinsic neuronal bursting has been observed, starting with the studies of Thoby-Brisson and Ramirez (2001) and Pena et al. (2004). Yes, our model by careful analysis of the sources of calcium flux activating ICAN, can explain the observations of Pena et al. of large perturbations of population activity amplitude with relatively smaller perturbations of oscillation frequency of the isolated pre-BötC in vitro, as well as similar results from the extensive experimental analysis of Koizumi et al. (2018) directly addressing this problem.

“The authors don't address the finding that blocking INaP alone has also no effect on the frequency (Pena et al., 2004), while others show that blocking INaP alone can block rhythmogenesis under control conditions. By contrast Pena et al. found that this is only the case in severe hypoxia. The modeling could explain these different experimental findings. For example the intracellular calcium concentrations or the network state could differ between different labs. Unfortunately, the authors don't pitch the paper like this, which seems like a missed opportunity. Thus, the authors should try to be less polarizing and report data as is (and as detailed below in more detail).”

Our statements are based largely on our published experimental results, including that application of riluzole in concentrations used by Pena et al. always led to cessation of the rhythm in rat slices under control conditions, which has also been reported in mouse slice as referenced (Toporikova et al., 2015). We have now noted in the Introduction “although some studies suggest that block of both INaP and ICAN are necessary to disrupt rhythmogenesis in vitro (Pena et al., 2004).” We have not attempted to solve this particular problem here, because of the focus in the modeling on tackling the problem of calcium dynamics and ICANactivation, but we have previously discussed potential differences in experimental preparations including network state that could possibly explain such differences in experimental results from a theoretical perspective in Jasinski et al. (2013) which we have now discussed (see Discussion section “Role of ICAN in the pre-BötC Respiratory Network”). We have also discussed how our model can explain the small perturbations of oscillation frequency accompanying large reductions in inspiratory population activity that have been observed experimentally with ICAN block by some investigators (e.g., Pena et al., 2004). We think that we have achieved a balanced and accurate presentation of experimental data throughout the present paper as the reviewers advise, and we have made revisions in response to the reviewers’ specific suggestions in this regard. The goal of the present study was to provide an explanation and formulate implications for the fact that the excitatory pre-BötC population is capable of varying the amplitude of its oscillations essentially independent of the frequency, which is analytically challenging, and our model formulation can potentially serve as a basis for resolving various experimental discrepancies with additional focused analyses.

“There is no need to try to push the agenda that only INaP is rhythmogenic when we know that ICAN has a big role in rhythmogenesis by regulating synaptic excitatory transmission.”

It is indeed the main conclusion of our study that ICAN contributes critically to excitatory synaptic interactions rather than intrinsic rhythmogenic properties of pre-BötC excitatory circuits. This is what explains why the network activity amplitude is very sensitive to manipulations of ICAN conductance while the oscillation frequency is not. We note that complete blockade of ICAN in our model can abolish synchronized network oscillations due to weakened excitatory synaptic transmission, which results in neuronal de-recruitment and desynchronization of the network, particularly when synaptic strength is low. So ICAN plays an important role in network activity synchronization but this is not rhythmogenesis per se. We have now noted this role of ICAN in synchronization in the Discussion section “Role of ICAN in the pre-BötC Respiratory Network”.

2) The authors should add a discussion regarding the neuron model and the robustness of the network activity. First, the model generates bursts of action potentials having an amplitude of around 30-40mV, as shown in Figure 4 and Figure 10. But in the experimental data of Koizumi et al., the amplitude of action potentials is more about roughly 80-100 mV (Figure 3 of the paper). Why is there such a big difference between the two? The spike height is indeed often used as a major criterion to assess the validity of a computational model at the cellular level. We agree that it is not a critical feature of the model used in this paper, but the issue deserves to be discussed. In particular, what happens to the bursting pattern if you attempt to reach a physiologically plausible spike height by increasing the transient sodium conductance (which might also require an increase in delayed-rectifier potassium conductance)? We suspect that bursting would be very fragile to an increase in Na conductance, due to the specific dynamical mechanisms at play in the model.

While it is certainly true that neuronal bursting capability could be affected by an increase in the transient sodium conductance without adjustments to other conductances, the bursting neuron model that we employ represents a stable bursting solution with the specific dynamical mechanisms at play. Indeed, complete biophysically accurate modeling of the respiratory neurons was not the goal of this study. Network synaptic interactions formulated in our model depend on times of spike occurrence and not on spike shape, height, or duration. Network activity is also characterized in terms of the number of generated spikes and not individual voltage dynamics. This allows us to claim that the results presented hold regardless of spiking characteristics. We, therefore, do not think that the discussion of the role of intrinsic spike properties would be especially important, which the reviewers seem to acknowledge. We also note that the transient sodium and delayed-rectifier potassium channel conductance (gKdr) dynamics in the current study are borrowed from previously published pre-BötC neuron models (Rybak et al., 2007; Smith et al., 2007; Jasinski et al., 2013). As in these previous studies, the spike height in the current model is approximately 45 mV which is smaller than the spikes shown in (Koizumi et al., 2018), although the spike height in the current model is within the range (35-80mV) of those shown in the numerous publications from other labs (e.g., Thoby-Brisson et al., 2001; Pena et al., 2004; Del Negro et al., 2005; Morgado-Valle et al., 2008; Ben-Mabrouk, 2010). We note that in the original model formulation of Butera et al. (1999a), INaP-dependent rhythmic bursting occurs with larger amplitude spikes than the present model, so bursting solutions exist with spike heights closer to those observed by Koizumi et al. We also note that in Figure 4E, F the spikes shown were deliberately truncated, as originally noted in the figure legend. Overall, we have attempted to demonstrate the robustness of the model network activity, particularly with respect to effects of changes in gCAN, in Figures 5 and 6. We have now indicated in the text the dependence of the rhythmogenic kernel on other factors such as leak conductance and gKdr in the section of the Discussion on “INaP-dependent rhythmogenic kernel” to give the readers some sense for what determines intrinsic rhythmic bursting in the model.

3) Along the same lines, Figure 4A suggests that the region of network bursting is quite narrow as compared to the region of network tonic spiking (especially if one increases the parameter range). What is the robustness of this bursting region against variability/heterogeneity in the other conductance values? Do you think that this result implies that physiological mechanisms need to tightly regulate the values of the conductances to maintain a viable network bursting activity (for instance a specific balance between synaptic strength and ICAN conductance), or could it be possible that this fragility is an artifact of the computational model?

We thank the reviewers for pointing out this issue. We added new text to the section (“Manipulating g-CAN in the CaSyn model is qualitatively equivalent to changing the strength of synaptic interactions”) describing Figure 4A to clarify the following. Because ICAN acts as a synaptic amplifier, the effective strength of synaptic interactions is roughly proportional to a product of the synaptic weight (SynW, shown as Synaptic Strength in Figure 4) andg-CAN. A transition from bursting to tonic spiking occurs when this effective excitation exceeds a certain critical value Wc in any imaginable network with recurrent excitation. That is why the bifurcation curve representing this boundary (between yellow and black) on the parameter plane (g-CAN, SynW) looks like a hyperbola g-CAN x SynW=Wc. Indeed, the stronger the basic synaptic interactions, the less room ICAN has to further facilitate them before the transition to the tonic regime occurs. So, all that Figure 4A implies is that the product of SynW and g-CAN should not exceed a certain value for the bursting to persist. It does not imply any specific balance between the two or non-robustness of the bursting regime.

4) The authors should succinctly describe the model at the beginning of the Results section. Indeed, it was difficult for us to understand the configuration of the model without having to go through the Materials and methods section, which is at the end of the manuscript and very detailed. It is important that the reader has a clear view of how neurons are modeled (ion channels), how many neurons constitute the circuit, how you compute "circuit activity amplitude", etc. before going through the analyses. Please also make sure to define every specific term when it is introduced. For instance, when describing how network probability is computed, Wmax is not defined, and it cannot be easily related to the "average weight of synaptic connections" without again going through the Materials and methods section. Also, you mention in the Abstract that ICAN conductance is distributed randomly, but not in the Materials and methods section. Which one is correct?

A succinct model description has now been added to the beginning of the Results sections and definitions for circuit activity amplitude and other terms such as synaptic strength have been added to the text when they are introduced. In the model, we also distributed ICAN conductance (g-CAN) based on a uniform distribution and this is now indicated in Table 1 and described in the Materials and methods. We actually found that the simulation results were not dependent on such a distribution as now stated in Materials and methods, and for the representative simulations presented in the figures, g-CAN values indicated in the figures represent the mean values for the simulated populations. Information about random ionic and synaptic conductance distributions has been removed from the Abstract, where such information is not really needed.

https://doi.org/10.7554/eLife.41555.018

Article and author information

Author details

  1. Ryan S Phillips

    1. Cellular and Systems Neurobiology Section, National Institute of Neurological Disorders and Stroke, National Institutes of Health, Bethesda, United States
    2. Department of Physics, University of New Hampshire, Durham, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8570-2348
  2. Tibin T John

    Cellular and Systems Neurobiology Section, National Institute of Neurological Disorders and Stroke, National Institutes of Health, Bethesda, United States
    Contribution
    Data curation, Formal analysis
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6825-7166
  3. Hidehiko Koizumi

    Cellular and Systems Neurobiology Section, National Institute of Neurological Disorders and Stroke, National Institutes of Health, Bethesda, United States
    Contribution
    Data curation, Formal analysis
    Competing interests
    No competing interests declared
  4. Yaroslav I Molkov

    1. Department of Mathematics and Statistics, Georgia State University, Atlanta, United States
    2. Neuroscience Institute, Georgia State University, Atlanta, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0862-1974
  5. Jeffrey C Smith

    Cellular and Systems Neurobiology Section, National Institute of Neurological Disorders and Stroke, National Institutes of Health, Bethesda, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Writing—original draft, Writing—review and editing
    For correspondence
    smithj2@ninds.nih.gov
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7676-4643

Funding

The Jayne Koskinas Ted Giovanis Foundation for Health and Policy

  • Ryan S Phillips

National Institute of Neurological Disorders and Stroke (Intramural Research Program of the National Institutes of Health)

  • Ryan S Phillips
  • Tibin T John
  • Hidehiko Koizumi
  • Jeffrey C Smith

National Institutes of Health (R01 AT008632)

  • Yaroslav I Molkov

National Institutes of Health (U01 EB021960)

  • Yaroslav I Molkov

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported in part by the Jayne Koskinas and Ted Giovanis Foundation for Health and Policy, the Intramural Research Program of the National Institutes of Health (NIH), National Institute of Neurological Disorders and Stroke (NINDS), and NIH Grants R01 AT008632 and U01 EB021960.

Senior Editor

  1. Ronald L Calabrese, Emory University, United States

Reviewing Editor

  1. Jan-Marino Ramirez, Seattle Children's Research Institute and University of Washington, United States

Reviewers

  1. Jan-Marino Ramirez, Seattle Children's Research Institute and University of Washington, United States
  2. Guillaume Drion

Publication history

  1. Received: September 4, 2018
  2. Accepted: February 7, 2019
  3. Version of Record published: March 25, 2019 (version 1)

Copyright

© 2019, Phillips et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 477
    Page views
  • 72
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

  1. Further reading

Further reading

    1. Computational and Systems Biology
    2. Developmental Biology
    Berta Verd et al.
    Research Article Updated
    1. Computational and Systems Biology
    2. Genetics and Genomics
    Dylan Kotliar et al.
    Tools and Resources Updated