The following features of the model outputs from new simulations (Figure 1) formed the basis for the design of experiments for model simulation-experimental data comparisons in the current study. Simulations were performed with the model specified in Phillips et al., 2019, with modifications as described in Materials and methods.

Figure 1 with 1 supplement see all

Download asset Open asset

Figure 1—source data 1

Related to Figure 1A–D.

https://cdn.elifesciences.org/articles/74762/elife-74762-fig1-data1-v1.xlsx

Download elife-74762-fig1-data1-v1.xlsx

1. The inspiratory rhythm in vitro is generated by a functionally distinct subpopulation of preBötC excitatory neurons with subthreshold-activating I_NaP that confers voltage-dependent rhythmic burst frequency control, providing a mechanism for frequency tuning by the regulation of baseline membrane potentials. Progressively depolarizing/hyperpolarizing neurons across the population by varying an applied current (I_App) increases/decreases population burst frequency over a wide dynamic range defined by the frequency tuning curve for the population (i.e. the relationship between applied current/baseline membrane potential and network bursting frequency).

2. Reducing neuronal persistent sodium conductance (g_NaP) decreases the population-level bursting frequency and alters the voltage-dependent rhythmogenic behavior with a shift in the frequency tuning curve such that the frequency dynamic range is reduced, consistent with previous simulations (Phillips and Rubin, 2019). Population amplitude is not as strongly reduced by decreases in g_NaP. Rhythm generation can be terminated at a sufficiently low level of g_NaP under baseline conditions of network excitability in vitro.

3. Reducing neuronal calcium-activated non-selective cation conductance [g_CAN] very strongly decreases population activity amplitude and, in contrast to reducing g_NaP, has little effect on population bursting frequency under baseline conditions, but strongly augments bursting frequency due to a shift in the frequency tuning curve that extends the frequency range at higher levels of population depolarization.

Model simulations showing the relations between I_App, population burst frequency, and burst amplitude at different levels of g_NaP and g_CAN are displayed in Figure 1, which illustrate these basic predictions of the model. Importantly, reducing the synaptic connection probability does not qualitatively change these predictions when the synaptic strength remains constant (see Figure 1—figure supplement 1), which allows for a fair assessment of the sensitivity of the network behavior to changes in connection probability (see also Phillips et al., 2019). The contrasting predictions for effects of g_NaP or g_CAN on network dynamics provide a clear basis for experimental testing of the model.

For comparisons with modeling results, we used electrophysiological recording from the preBötC inspiratory population and the hypoglossal (XII) motor output, in conjunction with pharmacological and optogenetic manipulations of preBötC glutamatergic neurons in the in vitro slices. Slices were prepared from transgenic mice that selectively express Channelrhodopsin-2 (ChR2) in glutamatergic neurons by Cre-dependent targeting via the vesicular glutamate transporter type-2 (VgluT2) promoter (Figure 2). ChR2-based photostimulation of preBötC glutamatergic neurons bilaterally enabled graded control of the baseline depolarization of this excitatory population to define the frequency tuning curves and associated population bursting amplitudes at different levels of pharmacological block of g_NaP or g_CAN.

Figure 2

Download asset Open asset

Figure 2—source data 1

Related to Figure 2D.

https://cdn.elifesciences.org/articles/74762/elife-74762-fig2-data1-v1.xlsx

Download elife-74762-fig2-data1-v1.xlsx

For optogenetic photostimulation, we established and validated a triple transgenic VgluT2-tdTomato-ChR2-enhanced yellow fluorescent protein (EYFP) mouse line, produced by crossing a previously validated VgluT2-Cre-tdTomato mouse (Koizumi et al., 2016) with a Cre-dependent ChR2-EYFP mouse strain. We verified Cre-driven ChR2-EYFP expression in VgluT2-positive tdTomato-labeled neurons within the preBötC region by two-photon microscopy in live neonatal mouse medullary slice (n=6) preparations in vitro, documenting heavy expression of the ChR2-EYFP fusion protein in processes and somal membranes of VgluT2-tdTomato-labeled neurons (Figure 2A). We also quantified by whole-cell patch-clamp recordings the photostimulation-induced depolarization of rhythmically active, ChR2-expressing inspiratory glutamatergic neurons (n=12 neurons from four slice preparations from VgluT2-tdTomato-ChR2- EYFP mouse line) at different levels of photostimulation (473 nm of variable laser power ranging from 0.5 to 5 mW) in the rhythmically active in vitro neonatal medullary slice preparations. These in vitro slices effectively isolate the bilateral preBötC along with hypoglossal (XII) motoneurons and nerves to monitor inspiratory XII activity (Koizumi et al., 2013; Koizumi et al., 2016), allowing laser illumination of the preBötC region unilaterally or bilaterally as well as whole-cell patch-clamp recordings from rhythmically active inspiratory preBötC neurons during laser illumination (Figure 2B). We performed targeted current-clamp recordings from rhythmic tdTomato-labeled VgluT2-positive neurons, in which co-expression of ChR2-EYFP fusion protein in the cell membranes was confirmed with two-photon, single optical plane live images (Figure 2B, right). Representative examples of ChR2-mediated membrane depolarization (Figure 2C) show that the membrane potential of a functionally identified VgluT2-positive preBötC inspiratory neuron was depolarized by ~6.0 mV at 1 mW, ~8.5 mV at 2 mW, and ~11.5 mV at 5 mW laser illumination. The light-induced depolarization had fast kinetics, occurring within ~50 ms and recovery within ~150 ms after terminating illumination. Summary data (n=12 neurons, mean ± SEM, Figure 2D) illustrates the laser-power-dependent ChR2-mediated depolarization of VgluT2-positive inspiratory preBötC neurons (4.73±0.24 mV at 0.5 mW, 6.01±0.21 mV at 1 mW, 8.55±0.29 mV at 2 mW, 10.21±0.34 mV at 3 mW, 11.13±0.3 mV at 4 mW, and 11.6±0.35 mV at 5 mW). The membrane depolarizations are significantly different (p<0.01 by ANOVA test comparing these groups). The subsequent Tukey’s HSD (honestly significant difference) test for pairwise comparison shows a significant difference with p<0.01 in all cases except p=0.4542 for 4 mW vs 5 mW. These results enabled us to define the relation between laser power and ChR2-mediated membrane depolarization for the inspiratory glutamatergic neuron population.

We performed site-specific bilateral photostimulations of the preBötC VgluT2-positive neuron population (Figure 3A) by sustained laser illumination (20 ms pulses at 20 Hz, 473 nm, 30–90-s pulse train durations) at graded laser powers between 0.25 and 2.0 mW (n=31 slice preparations). To systematically analyze relations between inspiratory burst frequency/burst amplitude and laser power, we applied single epochs of laser illumination and allowed for recovery from each epoch before changing the laser power (representative examples with simultaneous population activity recordings from the preBötC and XII are shown in Figure 3B, and summary data [mean ± SEM] is shown in Figure 3C). Bilateral photostimulation (0.25–2 mW) of the preBötC caused rapid and reversible increases of inspiratory burst frequency in a laser-power-dependent manner (85±7% increase at 0.25 mW, 174±7% at 0.5 mW, 262±6% at 1 mW, 326±6% at 2 mW, Wilcoxon matched-pairs signed rank test, p<0.01 in all cases). After termination of laser illumination, an inhibition of rhythmicity followed and subsequently recovered to the control frequency. The average recovery time (seconds) at each laser power was 48±4 at 0.25 mW, 88±4 at 0.5 mW, 128±9 at 1 mW, 182±11 at 2 mW, respectively. Concurrently, there was a significant laser power-dependent decrease of burst amplitude of inspiratory preBötC population and inspiratory XII activity during photostimulation (8±2% decrease at 0.25 mW, 15±2% at 0.5 mW, 25±3% at 1 mW, 44±4% at 2 mW, p<0.01 in all cases by Wilcoxon signed rank test). Bilateral photostimulation also induced tonic activity in preBötC population recordings, which was of higher amplitude at higher laser power, but did not consistently induce tonic activities in XII nerve output. Otherwise, XII burst amplitude and frequency always directionally mirrored preBötC glutamatergic population activity.

Figure 3

Download asset Open asset

Figure 3—source data 1

Related to Figure 3C.

https://cdn.elifesciences.org/articles/74762/elife-74762-fig3-data1-v1.xlsx

Download elife-74762-fig3-data1-v1.xlsx

These results demonstrate that preBötC glutamatergic neurons play an important role in regulating the inspiratory rhythm with a membrane voltage-dependent frequency control mechanism in the neonatal medullary slices in vitro. We note that in control experiments with rhythmic slice preparations from the VgluT2-tdTomato, non-ChR2-expressing transgenic mice (n=6), we tested for photo-induced perturbations of inspiratory burst frequency and amplitude, and confirmed that there were no significant perturbations of frequency and amplitude as a function of laser power (0.25–2 mW; ANOVA tests for frequency and burst amplitude, p=0.545 and 0.761, respectively).

After defining the voltage-dependent regulation of burst frequency and amplitude as described above, we investigated how pharmacological block of I_NaP in preBötC glutamatergic neurons affects this voltage-dependent control in the rhythmically active in vitro medullary slice preparations from the VgluT2-tdTomato-ChR2 mouse line. We first analyzed the pharmacological profile of block of I_NaP in preBötC inspiratory glutamatergic neurons with the I_NaP blockers tetrodotoxin citrate (TTX) at very low concentrations (5–20 nM, n=8) and also riluzole (5–20 µM, n=6) bath applied in the neonatal mouse medullary slice preparations as we have previously done for neonatal rat preBötC inspiratory neurons (Koizumi and Smith, 2008), but now analyzed for genetically identified mouse glutamatergic neurons. Whole-cell voltage-clamp recording from optically identified VgluT2-tdTomato-expressing inspiratory neurons was used to obtain current-voltage (I-V) relationships by applying slow voltage ramps (30 mV/s; –100 to +10 mV) and we computed TTX- and riluzole-sensitive I_NaP by subtracting I-V curves obtained before and after block of I_NaP with TTX or riluzole (Figure 4—figure supplement 1). Both TTX and riluzole attenuated the peak I_NaP amplitude (measured at ~–40 mV) in a dose-dependent manner. TTX reduced peak I_NaP by 36–63% at 5 nM, by 74–94% at 10 nM, and completely blocks I_NaP at 20 nM. Riluzole reduced the peak I_NaP by 36–66% at 5 µM, by 82–95% at 10 µM, and completely blocked I_NaP at 20 µM. In this dataset, there were five endogenous burster inspiratory neurons, which had significantly higher g_NaP/Cm (105.9±6.5 pS/pF) compared to non-endogenous burster inspiratory neurons (49.6±3.5 pS/pF) (p<0.05 by non-parametric Kolmogorov-Smirnov test), comparable to results from a previous study of neuronal properties performed in neonatal rat slices in vitro (Koizumi and Smith, 2008). However, the I_NaP attenuation by the blockers was not significantly different for these two groups (TTX: p=0.464 at 5 nM, p=0.982 at 10 nM, p=0.857 at 20 nM; riluzole: p=0.933 at 5 µM, p=0.4 at 10 µM, p=0.933 at 20 µM by non-parametric Kolmogorov-Smirnov tests). These pharmacological profiles are also comparable to the data obtained in neonatal rat preBötC inspiratory neurons in our previous in vitro studies (Koizumi and Smith, 2008). We note that these pharmacological profiles were obtained for neurons recorded at depths up to 50–75 µm in the slices; neurons located deeper in the slice tissue may exhibit less reduction of I_NaP at the time points when these more superficial neurons were analyzed due to tissue penetration-related, spatiotemporal non-uniformity of I_NaP attenuation.

We then comparatively analyzed control of XII inspiratory burst frequency and amplitude by bilateral preBötC photostimulation (473 nm, 20 Hz pulses, 30–90 s, trains of sustained photostimulation as described above) under I_NaP pharmacological blockade. As shown in Figure 4A and B, partial block of I_NaP (5–10 nM TTX, or 5–10 µM riluzole) gradually reduced XII inspiratory burst frequency, slightly reduced amplitude, and terminated inspiratory rhythmic bursting in almost all cases (n=6/7 at 5 nM TTX, n=8/8 at 10 nM TTX, n=6/7 at 5 µM riluzole, n=8/8 at 10 µM riluzole). In the cases where rhythmic bursting was terminated, after 30 min, we again bilaterally photostimulated the preBötC, which restored rhythmic bursting with voltage-dependent frequency control (0.25–2.0 mW) in all cases (n=14/14 with TTX, n=14/14 with riluzole). The summary data (n=6 slice preparations for TTX 5 nM, n=8 for TTX 10 nM, n=6 for riluzole 5 µM, n=8 for riluzole 10 µM, mean ± SEM plotted) in Figure 4C (left) shows that the inspiratory burst frequency significantly increases as a function of laser power (Wilcoxon signed rank test, p<0.05 in all cases). The frequency tuning curves (control, TTX 5 nM, TTX 10 nM, riluzole 5 µM, and riluzole 10 µM) are significantly different by ANOVA test (p<0.01). The tuning curves with TTX and riluzole shifted downward relative to control (TTX 5 nM vs control, TTX 10 nM vs control, riluzole 5 µM vs control, riluzole 10 µM vs control; in all cases, p<0.01 by subsequent Tukey’s HSD test), and more downwardly shifted with higher drug concentrations (Tukey’s HSD test, TTX 5 nM vs TTX 10 nM, p<0.01; riluzole 5 µM vs riluzole 10 µM, p<0.01). The frequency-laser power relations are not significantly different for the TTX and riluzole groups (Tukey’s HSD test, p=0.927 for TTX 5 nM vs riluzole 5 µM; p=0.969 for TTX 10 nM vs riluzole 10 µM).

Figure 4 with 1 supplement see all

Download asset Open asset

Figure 4—source data 1

Related to Figure 4C.

https://cdn.elifesciences.org/articles/74762/elife-74762-fig4-data1-v1.xlsx

Download elife-74762-fig4-data1-v1.xlsx

Summary data (n=6 slice preparations at TTX 5 nM, n=8 at TTX 10 nM, n=6 at riluzole 5 µM, n=8 at riluzole 10 µM, mean ± SEM plotted) of the relations between normalized inspiratory burst amplitude and laser power presented in Figure 4C (right) show a significant laser-power-dependent decrease of inspiratory burst amplitude with I_NaP attenuation (Wilcoxon signed rank test, p<0.05 in all cases). Although significant, changes in amplitude relative to control are smaller with TTX block of I_NaP. The amplitude vs laser power relation curves with partial I_NaP blockade are significantly different (p<0.01 by ANOVA). The curves with I_NaP attenuation (except for TTX 5 nM) are shifted significantly downward compared to control conditions (Tukey’s HSD test, TTX 5 nM vs control, p=0.304; TTX 10 nM vs control; riluzole 5 µM vs control, riluzole 10 µM vs control, p<0.05), but there were no significant differences for the relations with different concentrations of TTX or riluzole (TTX 5 nM vs TTX 10 nM, p=0.986; riluzole 5 µM vs riluzole 10 µM, p=0.259 by Tukey’s HSD test). The curves obtained under riluzole were more significantly downward-shifted than those under TTX (TTX 5 nM vs riluzole 5 µM, p<0.01; TTX 10 nM vs riluzole 10 µM, p<0.01 by Tukey’s HSD test).

Under conditions of complete block of I_NaP (>30 min after either 20 nM TTX or 20 µM riluzole), bilateral sustained preBötC photostimulation could not reinitiate the inspiratory rhythm at any laser power ranging from 0.25 to 5 mW (n=18/18 under TTX, n=11/11 under riluzole; Figures 4 and 5), but induced only graded tonic spiking in preBötC neurons with progressive increases of laser power >1 mW after block of I_NaP. This tonic activity without rhythmic bursting was also verified at the preBötC excitatory population activity level (Figure 5).

Figure 5

Download asset Open asset

Figure 5—source data 1

Related to Figure 5C and D.

https://cdn.elifesciences.org/articles/74762/elife-74762-fig5-data1-v1.xlsx

Download elife-74762-fig5-data1-v1.xlsx

We performed power spectrum analyses to determine if there are rhythmogenic mechanisms inherent within the preBötC excitatory network that do not depend on I_NaP and emerge with the graded levels of tonic excitatory population activity at the various levels of sustained bilateral photostimulation. These spectral analyses of the recorded graded tonic activity after block of I_NaP indicated that there was no rhythmic activity at each level of the photostimulation-induced tonic activity in comparison to the control activity. In control conditions there is clear definition of spectral peaks for the fundamental and higher harmonic frequencies as given by the Fast Fourier Transform of the rectified preBötC activity signal. This power spectrum reflects the periodicity and (non-sinusoidal) shape of the bursts of population activity, but the spectrum is flat after block of I_NaP (Figure 5, n=5), reflecting a complete absence of rhythmic activity. At the cellular level, current-clamp and voltage-clamp recordings confirmed the absence of neuronal rhythmic synaptic drive potentials/currents in preBötC inspiratory glutamatergic neurons (Figure 6).

Figure 6

Download asset Open asset

These results demonstrate that I_NaP in preBötC glutamatergic neurons is critically involved in generating the inspiratory rhythm with a voltage-dependent frequency control mechanism in vitro. The results also indicate that there are no apparent rhythmogenic mechanisms inherent within the excitatory network that can emerge at various levels of excitation after eliminating the I_NaP-dependent mechanisms.

We next investigated contributions of I_CAN/TRPM4 to the voltage-dependent control of XII inspiratory burst frequency and amplitude for comparisons to model predictions. We analyzed perturbations of burst frequency and amplitude by bilateral preBötC photostimulation (as described above) after pharmacologically blocking TRPM4 with a specific inhibitor of TRPM4 channels (9-phenanthrol, a putative blocker of I_CAN). Bath application of 9-phenanthrol at a concentration proposed to be selective for TRPM4 (50 µM) to the rhythmically active slice preparations from VgluT2-tdTomato-ChR2 neonatal mice (Figure 7A) significantly reduced the amplitude of inspiratory XII bursts (43±4% reduction, n=8, p<0.01 by Wilcoxon signed rank test), with little effects on burst frequency (1±1% increase, n=8, Wilcoxon signed rank test, p=0.461) in these slices, as we have previously described (Koizumi et al., 2018). Under these TRPM4 block conditions (steady-state at >60 min after 9-phenanthrol application, Figure 7C), photostimulation caused a significant increase of inspiratory burst frequency compared to control conditions at a given laser power (0.25–2.0 mW, Figure 7D black lines, Wilcoxon signed rank test, p<0.01 in all cases), accompanied by a significant laser-power-dependent decrease of burst amplitude compared to control (Figure 7D red lines, Wilcoxon signed rank test, p<0.01 in all cases). The relations between normalized inspiratory burst frequency and amplitude vs laser power summarized in Figure 7D (n=8 slice preparations, mean ± SEM plotted) show that TRPM4 block significantly shifts the frequency tuning curve upward compared to control conditions (black lines in Figure 7D, Wilcoxon signed rank test, p<0.01). The relation between burst amplitude and laser power is shifted downward significantly compared to control conditions (red lines in Figure 7D, Wilcoxon signed rank test, p<0.01). Note that under baseline conditions of excitability, TRPM4 block does not alter burst frequency, and the upward shift in the frequency tuning curve is opposite to the shift seen with TTX or riluzole block of I_NaP. These findings are consistent with the model predictions presented in Figure 1.

Figure 7

Download asset Open asset

Figure 7—source data 1

Related to Figure 7D.

https://cdn.elifesciences.org/articles/74762/elife-74762-fig7-data1-v1.xlsx

Download elife-74762-fig7-data1-v1.xlsx

For comparisons of model predictions and results from the experimental photostimulation and pharmacological manipulations (Figure 8), we incorporated in the original model (Phillips et al., 2019) a ChR2 current (${\mathrm{I}}_{\mathrm{C}\mathrm{h}\mathrm{R}2}$) to represent I_App for the excitatory neuron population and mimic effects of photo-induced neuronal depolarization. ChR2 was modeled by a four-state Markov channel (see Figure 9), based on biophysical representations of this channel in the literature (Williams et al., 2013), which we found could be parameterized as described in Materials and methods to closely match our experimental data on relations between laser power and membrane depolarization of preBötC inspiratory glutamatergic neurons (Figure 8A). To simulate TTX and riluzole block of I_NaP, we incorporated their distinct pharmacological mechanisms of action, as well as off-target effects of riluzole. TTX directly obstructs the Na⁺ pore (Hille, 2001), whereas riluzole shifts I_NaP inactivation dynamics in the hyperpolarizing direction (Hebert et al., 1994; Song et al., 1997; Ptak et al., 2005). Additionally, at the concentrations used for I_NaP block in the preBötC (≤20 µM), riluzole attenuates excitatory synaptic transmission (MacIver et al., 1996; Bellingham, 2011). Therefore, as in previous work (Phillips and Rubin, 2019), TTX application was simulated by a direct reduction in g_NaP, while a hyperpolarizing shift in V_h1/2 and a reduction in the conductance of recurrent excitatory synaptic connections by 20 and 25% were used to simulate 5 µM and 10 µM riluzole application, respectively. In agreement with the experimental results (Figure 4), the model simulations show that TTX or riluzole block of I_NaP decreases the population-level bursting frequency and amplitude at each level of laser power (Figure 8B–C). This reduction also alters the voltage-dependent rhythmogenic behavior, with a downward shift in the frequency tuning curve such that the frequency dynamic range is reduced and there is a complete termination of rhythm generation at a sufficient level of I_NaP block under baseline conditions of network excitability (i.e. laser power = 0 mW). Due to the simulated off-target attenuating effects of riluzole on excitatory synaptic transmission, the downward shift in the amplitude vs laser power curve is larger with riluzole than TTX application (as also shown in Phillips and Rubin, 2019) as seen in the experimental data.

Figure 8

Download asset Open asset

Figure 8—source data 1

Related to Figure 8A–C.

https://cdn.elifesciences.org/articles/74762/elife-74762-fig8-data1-v1.xlsx

Download elife-74762-fig8-data1-v1.xlsx

Figure 9

Download asset Open asset

I_CAN block (g_CAN reduction) more strongly decreases population activity amplitude in the model simulations than does I_NaP block and, in contrast to I_NaP block, has little effect on population bursting frequency under baseline conditions but strongly augments bursting frequency at higher levels of excitability. As such, I_CAN block results in an upward shift and steeper slope in the frequency tuning curve that extends the frequency range at the higher levels of population depolarization, consistent with the experimental results (Figure 7). Thus, the model predictions are qualitatively consistent with the experimental data for major features of excitatory preBötC network behavior.

As in the previous paper (Phillips et al., 2019), the preBötC model network is constructed with N=100 synaptically coupled excitatory neurons. Neurons are simulated with a single compartment described in the Hodgkin-Huxley formalism. The addition of $I_{ChR2}$ is a new feature of the current model compared to Phillips et al., 2019. In the updated model, the membrane potential $V_m$ for each neuron is given by the following current balance equation:

where $C_m$ is the membrane capacitance, $I_{Na}$, $I_K$, $I_{Leak}$, $I_{NaP}$, $I_{CAN}$, ${\displaystyle I_{Ca}}$, $I_{Syn}$, and $I_{ChR2}$ are ionic currents through sodium, potassium, leak, persistent sodium, calcium-activated non-selective cation, voltage-gated calcium, synaptic and ChR2 channels, respectively. For a full description of all currents other than $I_{ChR2}$, see Phillips et al., 2019. The description of $I_{ChR2}$ was adapted from Williams et al., 2013 and given by the following equation:

where $g_{ChR2}$ is the maximal conductance, $O_1$ and $O_2$ are light-sensitive gating variables, $\delta$ is the ratio of the two open-state conductances, $V_m$ is the membrane potential, and $E_{ChR2}$ is the reversal potential for ChR2. All parameters are given in Table 1. $G\left(V_m\right)$ is an empirically derived voltage-dependent function given by:

ChR2 is simulated as a four-state Markov channel with two open $O_1$ and $O_2$ and two closed states $C_1$ and $C_2$ (Figure 9). Transitions from $C_1$ to $O_1$ and from $C_2$ to $O_2$ are light-sensitive and controlled by the parameter $Irr$ which represents optical power with units of milliwatts (mW).

The rates of change between states are given by the following differential equations:

which are consistent with the condition that

The empirically estimated transition rate (ms^–1) parameters ${\omega }_1$, ${\displaystyle {\omega }_2}$, ${\gamma }_1$, ${\displaystyle {\gamma }_2}$, ${\displaystyle {\gamma }_3}$, ${\displaystyle {\gamma }_4}$, and ${\displaystyle {\gamma }_5}$ are defined as follows:

and

$F$ and $P$ represent the photon flux (number of photons per molecule per second) and time- and irradiance-dependent activation function for ChR2, respectively. These quantities take the form

where ${\sigma }_{ret}$ is the absorption cross-section, $\lambda$ is the wavelength of max absorption, ${\omega }_{loss}$ is a scaling factor for loss of photons due to scattering/absorption, $h$ is the Planks constant, $C$ is the speed of light, and ${\tau }_{ChR2}$ is the time constant of ChR2 activation. Model parameter values are given in Table 1.

Initial model predictions presented in Figure 1 were generated without modification of the initial model tuning presented in Phillips et al., 2019. To match experimental data (Figure 8), parameters were slightly modified from the original model. The updated list of parameters is given in Table 1, where parameter adjustments are indicated with red font.

All simulations were performed locally on an eight-core computer running the Ubuntu 20.04 operating system. Simulation software was custom written in C++ and compiled with g++ version 9.3.0 (source code is provided as Source code 1). Numerical integration was performed using the exponential Euler method with a fixed step-size (dt) of 0.025 ms. In all simulations, the first 50 s of simulation time was discarded to allow for the decay of any initial condition-dependent transients.

The assumptions and neurobiological simplifications of the original model were previously elaborated in Phillips et al., 2019. We have extended the model to allow comparisons of our new experimental data and model simulation results. The updates included incorporating a ${\mathrm{I}}_{\mathrm{C}\mathrm{h}\mathrm{R}2}$ to represent I_App for the entire excitatory neuron population and to mimic effects of photo-induced neuronal depolarization. We modeled ChR2 by a four-state Markov channel based on channel biophysical representations in the literature (Williams et al., 2013). We found that this biophysical model could be parameterized to closely match our experimental data on relations between laser power and membrane depolarization of identified preBötC inspiratory glutamatergic neurons. Furthermore, we found that this channel model parameterized from our cellular-level data yielded frequency tuning curves for the bilateral excitatory population under baseline conditions and with pharmacological perturbations that were in close agreement with the experimentally measured tuning curves. We note that in the model, the assumption is that all excitatory network neurons have the same level of depolarization at a given laser power. In the experimental in vitro slice preparations, this may not occur if the neuronal expression levels of ChR2 are non-uniform. We have used transgenic strategies in mice to drive neuronal expression of ChR2 throughout the population of glutamatergic neurons, which has the advantage that there should be efficiency in coverage of this population in terms of channel expression. In support of this assumption, our single neuron electrophysiological data show a small SEM in the relations between neuronal depolarization and laser power, but it is not possible to fully quantify the uniformity of ChR2 expression levels across the entire population of glutamatergic inspiratory neurons.

We also note that the transgenic strategy used here may result in expression of ChR2 in glutamatergic terminals, including on preBötC inhibitory neurons that are proposed to interact with the excitatory neurons (e.g. Ausborn et al., 2018; Ashhad and Feldman, 2020), as well as fibers of passage from non-preBötC neurons, which could result in off-target effects. However, as discussed above, we have calibrated the ChR2 model to match the levels of depolarization at a given laser power that we have measured experimentally at the neuronal level. These measurements represent the total amount of neuronal depolarization which would include any contributions from photostimulated glutamatergic fibers and terminals. Thus, our ChR2 model and simulated photostimulation experiments adequately account for all sources contributing to neuronal depolarization and such off-target effects related to ChR2 expression, if present, would not impact the conclusions of this study.

To compare the model behavior with data from the pharmacological manipulations, we simulated blockade of I_NaP by TTX or I_CAN inhibition by 9-phenanthrol by reducing the channel conductances g_NaP and g_CAN, respectively. Blockade of I_NaP by riluzole was simulated by a hyperpolarizing shift in the inactivation parameter V_h1/2 and by a partial reduction in W_max (–20% and –25% for simulation of 5 µM and 10 µM riluzole), as described and justified previously (Phillips and Rubin, 2019) to take into account known or proposed pharmacological actions of riluzole from the literature. With this approach, the simulation results were directionally consistent with the data showing that I_NaP attenuation, by either low concentrations of TTX or by riluzole, and I_CAN blockade have opposite effects on the relationship between network excitability controlled by photostimulation and preBötC population burst frequency. We verified the presence of I_NaP in preBötC glutamatergic inspiratory neurons from voltage-clamp measurements and demonstrated partial and complete block of this current by TTX or riluzole under our in vitro experimental conditions. A difference between the modeling and experimental results that may represent a limitation of the model is that the model network is more sensitive to I_NaP block in terms of rhythm perturbation than suggested by the experimental results. In both model simulations and experiments, partial or complete block of this current can terminate rhythm generation; however, rhythm generation can stop at a smaller percentage reduction in I_NaP in the model, where it is assumed that I_NaP in all network neurons is uniformly reduced by the same amount at a given simulated percent block. This uniformity may not be the case in the in vitro network, where heterogeneity of cellular biophysical properties and spatial variability of pharmacological actions associated with drug penetration problems in the slice could cause non-uniformity of block at a given drug concentration.

Furthermore, the biophysical and pharmacological properties of I_NaP in preBötC inspiratory neurons are not fully characterized. Our recent biophysical characterization of I_NaP for these neurons has confirmed that this current is slowly inactivating (Yamanishi et al., 2018), which is a basic assumption of the present and previous models (Phillips et al., 2019), although the details of the voltage-dependent kinetic properties are more complex than represented by the first-order kinetics in the models. Moreover, some additional currents or dynamical processes (e.g. inhibitory currents from local circuit connections, synaptic depression, burst terminating currents, ionic dynamics; Thoby-Brisson et al., 2000; Hayes et al., 2008; Jasinski et al., 2013; Kottick and Del Negro, 2015; Phillips et al., 2018; Baertsch et al., 2018; Juárez-Vidales et al., 2021; Abdulla et al., 2022; Revill et al., 2021), which are not considered in the model, may play a role in augmenting/shaping inspiratory bursting, contribute to robustness of the rhythm such as the electrogenic transmembrane Na⁺/K⁺-ATPase pump (Jasinski et al., 2013), as well as affect preBötC dynamics on longer time scales than those investigated here. For example, the current model does not predict the minutes long decrease in network burst frequency that may occur following optogenetic stimulation in some experiments (see Figure 3, but also see Figure 7B, where this effect is less strong). Similar activity-dependent inhibitory effects in mammalian motor networks and tadpole locomotion neurocircuits have been shown to be mediated by the Na⁺/K⁺-ATPase pump (Picton et al., 2017; Picton et al., 2018). Therefore, capturing this feature in our model may require inclusion of the Na⁺/K⁺-ATPase pump which is beyond the scope of this study and, along with the exploration of other possible mechanisms that may impact or shape preBötC inspiratory bursting, is left for future investigations. Nonetheless, with the kinetic and pharmacological properties incorporated in our model, the directional trends in the experimental results for the various pharmacological manipulations of I_NaP and I_CAN are correctly predicted by the model simulations with very minor additional adjustments of key parameters relative to the initial model developed in Phillips et al., 2019.

Many previous experimental and theoretical analyses have focused on the rhythmogenic mechanisms operating under in vitro conditions in rhythmically active slices to provide insight into the biophysical and circuit processes involved, which as discussed in a number of reviews, has potential relevance for rhythm generation in vivo (Feldman and Del Negro, 2006; Richter and Smith, 2014; Phillips et al., 2019). There is an agreement that preBötC circuits have intrinsic autorhythmic properties, particularly because when isolated in slices in vitro, these circuits continue under appropriate conditions of excitability to generate rhythmic activity that drives behaviorally relevant inspiratory hypoglossal motoneuronal output. However, there is currently no consensus about the underlying rhythmogenic mechanisms. The endogenous rhythmic activity in vitro has been suggested to arise from various cellular and circuit biophysical mechanisms, including from a subset of intrinsically bursting neurons which, through excitatory synaptic interactions, recruit and synchronize neurons within the network (pacemaker-network models; Rekling and Feldman, 1998; Ramirez et al., 2004; Toporikova and Butera, 2011), or as an emergent network property involving recurrent excitation (Jasinski et al., 2013) and/or synaptic depression (group pacemaker model; Rubin et al., 2009a). More recent ‘burstlet’ models for rhythm generation emphasize how rhythm emerges in low excitability states due to synchronization of subsets of excitatory bursting neurons (Kam et al., 2013; Del Negro et al., 2018). Extensions of our current model (Phillips and Rubin, 2022) with our proposed I_NaP and I_CAN biophysical mechanisms for rhythm and activity amplitude generation can explain the generation of burstlets and synchronized population activity.

The present experimental results are consistent with previous findings indicating that the subthreshold-activating, slowly inactivating I_NaP is a critical neuronal conductance for inspiratory rhythm generation and neuronal bursting in vitro (Koizumi and Smith, 2008; Yamanishi et al., 2018). As our previous experimental and modeling results have indicated, the voltage-dependent and kinetic properties of these channels are fully capable of orchestrating rhythmic bursting at neuronal (e.g. Yamanishi et al., 2018) and excitatory neuron population levels (Koizumi et al., 2016). In the experiments blocking I_NaP, we tested by power spectral analyses for rhythmogenic mechanisms that do not depend on I_NaP after fully blocking this conductance as verified by our cellular-level measurements. These experiments were designed to determine if there are additional rhythmogenic mechanisms inherent within the excitatory network that can emerge during graded levels of tonic network excitation as controlled in our experiments by sustained bilateral photostimulation. We did not find any coherent rhythmic population activity under these conditions that would indicate the presence of other important rhythmogenic mechanisms capable of producing rhythmic population bursting on any time scale, which is consistent with our model. These results do not support the concept that rhythm generation under in vitro conditions is an emergent network property through recurrent excitation (Jasinski et al., 2013) and/or synaptic depression (group pacemaker model; Rubin et al., 2009a). Such emergent rhythms are theoretically possible as shown by previous modeling studies (Rubin et al., 2009a; Jasinski et al., 2013; Guerrier et al., 2015), but the regenerative population burst-generating and burst-terminating mechanisms incorporated in these models are apparently not sufficiently expressed in the in vitro preBötC network.

We note that our new measurements indicate that while I_CAN /TRPM4 activation is not involved in rhythm generation under baseline conditions of network excitability, since inhibiting this current does not significantly affect the rhythm, we also show that reducing this current does shift the frequency tuning curve at higher levels of network excitation. This occurs because activation of this current augments excitatory synaptic interactions, which can affect network bursting frequency. As discussed above, these observed effects of reducing I_CAN /TRPM4 on bursting frequency at elevated levels of network excitation, which we were able to reveal with our photostimulation paradigm, is a basic feature and entirely consistent with predictions of our model.

Our new experimental results confirm our previous experimental and modeling results indicating that endogenous activation of I_CAN /TRPM4 is critically involved in generating the amplitude of population activity but not the rhythm, which is also consistent with results from a genetic knockdown study of TRPM4 channels in preBötC neurons in vivo (Picardo et al., 2019). The correspondence between the experimental data and model predictions supports the concept in the model that activation of I_CAN is largely due to synaptically activated sources of neuronal Ca²⁺ flux such that I_CAN contributes to the excitatory inspiratory drive potential and regulates inspiratory burst amplitude by augmenting the excitatory synaptic current. Our data also show that I_NaP is involved, to a smaller degree, in generating the amplitude of rhythmic population activity at a given level of excitatory network excitation. This is due to the subthreshold activation of I_NaP and its voltage-dependent amplification of synaptic drive; indeed, application of riluzole, which is thought to impact synaptic transmission, affected population amplitude much more than application of TTX (Figures 4 and 8). This contribution to amplitude is smaller than that from I_CAN /TRPM4 activation, however. In general, our results support the concept from our original model that I_CAN activation in a subpopulation of preBötC excitatory neurons is critically involved in amplifying synaptic drive from a subset of neurons whose rhythmic bursting is critically dependent on I_NaP . This latter subpopulation forms the kernel for rhythm generation in vitro.

The role of I_NaP in rhythm generation within preBötC circuits is highly debated in the field with diverse and contradictory experimental results (Del Negro et al., 2002; Peña et al., 2004; Ramirez et al., 2004; Del Negro et al., 2005; Feldman and Del Negro, 2006; Smith et al., 2007; Pace et al., 2007; Koizumi and Smith, 2008; Ashhad and Feldman, 2020). Pharmacological studies (Smith et al., 2007; Koizumi and Smith, 2008) using riluzole or low concentrations of TTX to block I_NaP in preBötC neurons/circuits, as in the present studies, demonstrated a large reduction of preBötC inspiratory bursting frequency at cellular and network levels, with relatively smaller reductions of inspiratory network activity amplitude; fully blocking I_NaP completely eliminated inspiratory rhythm generation within isolated preBötC circuits in vitro. These experimental results are consistent with the computational hypothesis presented here and earlier that I_NaP is an essential biophysical component of inspiratory rhythm generation in vitro (Butera et al., 1999a; Butera et al., 1999b; Smith et al., 2000; Phillips et al., 2019; Phillips and Rubin, 2019).

However, this hypothesis has been challenged by some previous studies (Del Negro et al., 2002; Peña et al., 2004; Pace et al., 2007). In particular, Pace et al., 2007 found that microinjecting 10 µM riluzole or 20 nM TTX directly into the raphe obscurus but not the preBötC stopped inspiratory rhythm generation in vitro. Moreover, preBötC rhythm generation could be restarted by bath application of 0.5 µM substance P following bath application and microinjection of riluzole. These results suggest that bath application of riluzole or TTX does not compromise the fundamental mechanism(s) of preBötC rhythmogenesis but rather the level of neuronal excitability via off-target effects in the raphe obscurus. These results have not been reproduced, however, and in a follow-up study, Koizumi and Smith, 2008 observed that bath application or bilateral microinfusion in the preBötC of 10 µM riluzole or 20 nM TTX reliably abolished preBötC rhythm generation. The same is true for bath application of TTX in the in vitro preBötC island preparations that do not contain the raphe obscurus. Moreover, bath application of 1 µM substance P could only restart the preBötC rhythm when I_NaP block was incomplete. The present study is also a direct test of the conclusions reached by Pace et al., 2007. If, as suggested by Pace et al., 2007, bath application of TTX or riluzole impacts the inspiratory rhythm by reducing preBötC excitability rather than by affecting the essential mechanism(s) of rhythm generation, then increasing preBötC excitability via optogenetic stimulation should restart the rhythm even after complete I_NaP blockade. Our results show that, to the contrary, the preBötC is incapable of generating rhythmic output after complete I_NaP block even under optogenetic stimulation (Figures 4 and 5), demonstrating that I_NaP is essential for preBötC rhythm generation in this reduced in vitro preparation.

The off-target effects of riluzole are also a concern when interpreting its impact on preBötC rhythm and pattern generation. Specifically, at the concentrations used in these experiments (≤20 µM) riluzole is reported to attenuate excitatory synaptic transmission (MacIver et al., 1996; Bellingham, 2011). Reducing the strength of excitatory synaptic transmission has previously been shown to decrease preBötC network burst amplitude and slightly increase frequency in computational models (Phillips et al., 2019; Phillips and Rubin, 2019) and experimental studies (Johnson et al., 1994). Moreover, Phillips et al., 2019 showed that a 7–30% reduction in the strength of excitatory synaptic transmission was required to match the dose-dependent reduction of network burst amplitude seen with progressive riluzole infusion into the preBötC in the in vitro slice preparations. In the current study, the larger downward shift in the network amplitude vs laser power curve seen in the model with simulated riluzole application (compared to TTX), as also seen in our in vitro slice preparations (Figure 8B–C), is due to the 20% (at 5 µM) to 25% (at 10 µM) attenuation of excitatory synaptic conductances. Therefore, the subtly different effects of riluzole and TTX application on network burst amplitude in the current study support the idea that riluzole diminishes excitatory synaptic transmission by approximately these factors.

If I_NaP is critical for rhythm generation in the preBötC, then what can explain the findings (Del Negro et al., 2002; Peña et al., 2004; Pace et al., 2007) that seem to contradict this conclusion? Although more experimental work is needed to provide a definitive answer, possible explanations can be derived from computational modeling by considering non-uniformity of I_NaP blockade, in vitro slice thickness, and the order in which neurons (patterning vs rhythmogenic neurons) in the network are affected (Phillips and Rubin, 2019). Penetration of bath applied or microinjected I_NaP blockers in the in vitro slice preparations depends on passive diffusion. Therefore, the magnitude and progression of I_NaP block across the slice will not be uniform and drug penetration may be incomplete in thicker slices, even with microinjection within the preBötC; limited drug diffusion could explain why microinjection of TTX or riluzole into the preBötC may fail to stop rhythm generation (Pace et al., 2007). Moreover, the preBötC network dynamically expands/contracts in the rostralcaudal axis (Baertsch et al., 2019), which suggests that rhythmogenic neurons are preferentially located near the center of a larger population of pattern forming inspiratory neurons. This observation leads to the prediction that I_NaP block in thick slices via bath application of TTX or riluzole will primarily affect the amplitude of the inspiratory rhythm rather than its frequency and persistence, which matches the findings of in vitro experiments using thick slices (Del Negro et al., 2002; Peña et al., 2004; Pace et al., 2007). Another factor that impacts dynamics in thicker slices than those used in our studies is that they contain not only the preBötC, but also other network components that may introduce additional mechanisms. Indeed, with more intact brainstem respiratory pattern generation networks that exhibit more aspects of the full inspiratory and expiratory neuronal activity pattern in mature rodents, I_NaP blockers do not disrupt rhythm generation (Smith et al., 2007; Rubin et al., 2009b; Rubin and Smith, 2019). In the more intact system, neuronal dynamics in the preBötC are controlled by more complex synaptic interactions, including inhibitory circuit interactions, with new rhythmogenic regimes (Rubin et al., 2009b; Rubin and Smith, 2019). Interestingly, modeling results suggest that there are network states and pharmacological conditions under which the inspiratory rhythm may be highly dependent upon I_NaP in the more intact respiratory pattern generation network (Phillips and Rubin, 2019; Rubin and Smith, 2019), and these predictions require further experimental testing.

Regardless of this greater complexity, our present study has confirmed our previous results that I_NaP is essential for inspiratory rhythm generation in the isolated preBötC excitatory circuits in vitro, and represents a major rhythmogenic mechanism in this reduced state of the respiratory network. An important prediction of our model (see Phillips et al., 2019, Figure 9) is that the subpopulation of excitatory neurons with I_NaP-dependent rhythmogenic properties forming the rhythmogenic kernel may be relatively small compared to the neuronal subpopulation(s) with the I_CAN-dependent mechanism that critically amplifies the rhythmic drive from the kernel to generate and control the amplitude of excitatory population activity. This prediction also requires direct experimental testing, such as by dynamic multicellular Ca²⁺ imaging in glutamatergic neurons in vitro (Koizumi et al., 2018).

Our exploitation of optogenetic control of glutamatergic neuron network excitability, in combination with specific pharmacological manipulations of neuronal conductances, has enabled rigorous testing of predictions of our previous model of preBötC excitatory circuits. The basic predictions of the model for cellular and network behavior under the experimental conditions tested show a strong overall agreement with our experimental results. This agreement advances our understanding of neuronal and circuit biophysical mechanisms generating the rhythm and amplitude of inspiratory activity in the brainstem preBötC inspiratory oscillator in vitro and demonstrates the predictive power of our model.

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

This work was supported in part by the Intramural Research Program of the National Institutes of Health (NIH), National Institute of Neurological Disorders and Stroke (JS, HK), and grants from NSF DMS 1951095 (JR), NSF DMS 1724240 (JR, RSP), NIH R01 AT008632 (YM), NIH U01 EB021960 (YM) and Georgia State University B&B Seed Grant (YM).

This study was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health, the National Institute of Neurological Disorders and Stroke (Animal Study Proposal number: 1154-21).

This is an open-access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.