Predictions and experimental tests of a new biophysical model of the mammalian respiratory oscillator

Previously our computational modeling studies (Phillips et al., 2019) proposed that neuronal persistent sodium current (INaP) and calcium-activated non-selective cation current (ICAN) are key biophysical factors that, respectively, generate inspiratory rhythm and burst pattern in the mammalian preBötzinger complex (preBötC) respiratory oscillator isolated in vitro. Here, we experimentally tested and confirmed three predictions of the model from new simulations concerning the roles of INaP and ICAN: (1) INaP and ICAN blockade have opposite effects on the relationship between network excitability and preBötC rhythmic activity; (2) INaP is essential for preBötC rhythmogenesis; and (3) ICAN is essential for generating the amplitude of rhythmic output but not rhythm generation. These predictions were confirmed via optogenetic manipulations of preBötC network excitability during graded INaP or ICAN blockade by pharmacological manipulations in slices in vitro containing the rhythmically active preBötC from the medulla oblongata of neonatal mice. Our results support and advance the hypothesis that INaP and ICAN mechanistically underlie rhythm and inspiratory burst pattern generation, respectively, in the isolated preBötC.


Introduction
The cellular and circuit biophysical mechanisms generating the breathing rhythm critical for life in mammals have been under investigation for decades without clear resolution (Richter and Smith, 2014;Del Negro et al., 2018;Ramirez and Baertsch, 2018). The roles of various transmembrane ionic currents in this process have received considerable attention, including our recent updated computational modeling  of the biophysical mechanisms operating in the mammalian respiratory oscillator within the brainstem preBötzinger complex (preBötC), which contains essential rhythmogenic neuronal circuits.
Simulations from our updated model of preBötC excitatory circuits support distinct functional roles for cellular-level Na + -and Ca 2+ -based biophysical mechanisms that have been of intense experimental and theoretical interest for understanding the operation of preBötC circuits . These mechanisms rely on a slowly inactivating neuronal persistent sodium current (I NaP ) and a calcium-activated non-selective cation current (I CAN ) mediated by transient receptor potential M4 (TRPM4) channels coupled to intracellular calcium dynamics. The model proposes how these cellularlevel mechanisms operating in excitatory circuits can underlie two critical functions of preBötC circuits: generation of the inspiratory rhythm and regulation of the amplitude of inspiratory population activity. In essence, the model advances the concepts that (1) a subset of excitatory circuit neurons, whose rhythmic bursting in the isolated preBötC in vitro is critically dependent on I NaP , forms an excitatory neuronal kernel for rhythm generation, and (2) excitatory synaptic drive from the rhythmogenic kernel population is critically amplified by I CAN activation in a recruited subset of preBötC excitatory neurons to generate overall population activity of sufficient amplitude to induce downstream excitatory circuits to produce inspiratory motor output.
While these concepts have some experimental support (see discussions and model-experimental data comparisons in , the mechanisms incorporated in the model and their predicted effects on circuit behavior require more extensive experimental testing. Indeed, the basic rhythmogenic I NaP -based biophysical mechanism represented by the model remains controversial . The concept that activation of I CAN is largely due to synaptically activated sources of neuronal Ca 2+ flux, such that I CAN tunes inspiratory burst amplitude in proportion to the level of excitatory synaptic current from the kernel, also requires additional experimental support. In the present modeling and experimental study, to build on our previous work , we have performed experimental tests on the rhythmically active in vitro network in slices from neonatal transgenic mice using a combination of electrophysiological analyses, pharmacological perturbations of I NaP or I CAN , and optogenetic manipulations of the preBötC excitatory (glutamatergic) population involved. Our new comparisons of model simulations and experimental data provide further evidence for the predictive power and validity of the model, while also indicating some limitations.

Model concepts and predictions
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 , with modifications as described in Materials and methods.
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  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 . The contrasting predictions for effects of g NaP or g CAN on network dynamics provide a clear basis for experimental testing of the model.

Experimental design and results
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 .

Cre-dependent targeting of preBötC glutamatergic neurons for photostimulation experiments
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 frequency and amplitude of model neuronal population activity, and in (D) at fixed levels of reduction (10%, 20%) in g NaP from control. Color scale bars for values of burst frequency and amplitude are at right of plots in (A) and (C). Frequencies and amplitudes of population activity in (C) and (D) are normalized to control values.
The online version of this article includes the following source data and figure supplement(s) for figure 1: Source data 1. Related to Figure 1A-D.
Figure supplement 1. Comparison of model simulation predictions with a lower synaptic connection probability of 13% (p Syn =0.13), as experimentally approximated by Rekling et al., 2000, with those presented in Figure 1.

Perturbations of inspiratory rhythm and burst amplitude by photostimulation of glutamatergic neurons within the preBötC region
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. 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 photoinduced 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).
VgluT2-tdTomato neurons, as seen in the merged image. Abbreviations: d, dorsal; l, lateral; NAsc, semi-compact subdivision of nucleus ambiguus. (B) Overview of experimental in vitro rhythmic slice preparation from a neonatal VgluT2-tdTomato-ChR2 transgenic mouse showing whole-cell patchclamp recording from a functionally identified preBötC inspiratory VgluT2-positive neuron with unilateral preBötC laser illumination (0.5-5 mW) to test for neuronal expression of ChR2, and suction-electrode extracellular recordings from hypoglossal (XII) nerves to monitor inspiratory activity. Regional photostimulation was accomplished with a 100-µm diameter optical cannula positioned at the slice surface above the preBötC region. Abbreviations: V4, fourth ventricle; IO, inferior olivary nucleus. Upper right in (B) two-photon single optical plane images showing an imaged preBötC inspiratory neuron (arrow) targeted for whole-cell recording. From left, Dodt gradient contrast structural image, VgluT2-Cre driven tdTomato labeling, ChR2-EYFP expression, and merged image that confirm co-expression of tdTomato and ChR2-EYFP. The lower right traces are current-clamp recordings from the VgluT2-positive preBötC neuron shown in the above images, illustrating inspiratory spikes/bursts synchronized with integrated inspiratory XII nerve activity (∫XII). The traces shown at right within the dashed box are expanded time scale traces of the bursting activity within the dashed box in the left traces. (C) The membrane potential (Vm) of the neuron shown in (B) was depolarized during photostimulation by ~6.0 mV at 1 mW, ~8.5 mV at 2 mW, and ~11.5 mV at 5 mW of laser power (spikes are truncated). The neuron was hyperpolarized from resting baseline potential to -74 mV by applying constant current in this example to reveal the magnitude of the light-induced membrane depolarization. Photostimulation was performed in the interval between inspiratory population bursts. (D) Summary data (n=12 neurons from four slice preparations, mean ± SEM) showing the laser-power-dependent, ChR2-mediated neuronal membrane depolarization from baseline (ΔV M ) of VgluT2-positive preBötC inspiratory neurons.
The online version of this article includes the following source data for figure 2: Source data 1. Related to Figure 2D.  . Photostimulation of the bilateral preBötzinger complex (preBötC) vesicular glutamate transporter type-2 (VgluT2)-positive neuron population caused laser-power-dependent increases of inspiratory burst frequency and decreases of burst amplitude. (A) Overview of experimental in vitro rhythmically active slice preparation from neonatal VgluT2-tdTomato-channelrhodopsin-2 transgenic mouse with macro-patch electrodes on the preBötC region for recording preBötC population activity and suction electrodes on hypoglossal (XII) nerves to monitor inspiratory motor activity during bilateral preBötC laser illumination (0.25-2 mW, 473 nm wavelength) with optical cannula (100 µm diameter) positioned obliquely to illuminate the preBötC region. Abbreviations: NAsc, nucleus ambiguus semi-compact subdivision; V4, fourth ventricle; IO, inferior olivary nucleus. (B) Representative examples of epochs of bilateral preBötC laser illumination and effects on inspiratory burst frequency and burst amplitude. The upper traces show the integrated macro-patch recordings from the preBötC inspiratory neuron population (∫preBötC), middle traces show the integrated inspiratory XII activity (∫XII), and the bottom traces show the inspiratory burst frequency (time-based moving median in a 10 s window). The sustained laser illumination with the laser intensity is given by blue shading. Low-intensity illumination (0.5 mW) caused significant increase (~149% in this example) of inspiratory burst frequency and decrease of inspiratory burst amplitude of both preBötC and XII population activity (~22 and~21%, respectively, left panel). Higher intensity laser illumination (2 mW; right panel) caused a larger increase (~295%) of inspiratory frequency and decrease of burst amplitude of inspiratory preBötC and XII activity (~65 and~44%, respectively). Also note that photostimulation induced tonic activity (indicated by baseline shift) in ∫preBötC population activity recordings. (C) Summary data of relations between inspiratory burst frequency and amplitude vs photostimulation laser power (n=31 slices; data points plotted are mean ± SEM).
The online version of this article includes the following source data for figure 3: Source data 1. Related to Figure 3C.
Regulation of burst frequency and amplitude by I NaP in preBötC glutamatergic neurons analyzed with combined optogenetic and pharmacological manipulations 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.
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 Figure 4. Perturbations of inspiratory burst frequency and amplitude by bilateral preBötzinger complex (preBötC) photostimulation during pharmacological block of neuronal persistent sodium current (I NaP ). (A) Example recordings of integrated XII activity (∫XII) with bath application of low concentration of tetrodotoxin citrate (TTX; 5 nM), which gradually decreased inspiratory burst frequency and completely stopped the rhythm within ~25-30 min. The amplitude of XII activity also gradually decreased (~36% before the rhythm stopped in this example). (B) Under these conditions of partial block of I NaP , bilateral preBötC photostimulation (30 min after the rhythm stopped, blue-shaded epoch) could reinitiate the rhythm, which was also laser-power-dependent (~202% increase at 1 mW compared to the control, and ~306% increase at 2 mW in the example shown). (C) Summary plots (TTX 5 nM, n=6; TTX 10 nM, n=8; TTX 20 nM, n=18; riluzole [RZ] 5 µM, n=6; RZ 10 µM, n=8; RZ 20 µM, n=8; mean ± SEM plotted) of the relations between laser power and normalized inspiratory burst frequency indicate laser-power-dependent, significant increases of frequency in all cases. The curves for higher concentration of TTX and riluzole are downward-shifted compared to those for the lower concentration as well as those under control conditions (before drug applications). Right panel shows summary plots (TTX 5 nM, n=6; TTX 10 nM, n=8; TTX 20 nM, n=18; RZ 5 µM, n=6; RZ 10 µM; n=8; RZ 20 µM, n=8; mean ± SEM plotted) of the relations between laser power and normalized XII burst amplitude indicating laser-power-dependent, significant decreases of burst amplitude in all cases. The curves under the I NaP partial block (except for at TTX 5 nM) are downward-shifted compared to those under the control conditions. The decrease in burst amplitudes under RZ is more significant than those under TTX, although there are no significant differences between different concentrations of TTX or RZ. The loss of network rhythmic bursting activity after complete block of I NaP (20 nM TTX or 20 µM RZ) is reflected by the zero frequency (left) and amplitude (right) points for 20 µM RZ and 20 nM TTX on the plots.
The online version of this article includes the following source data and figure supplement(s) for figure 4: Source data 1. Related to Figure 4C.  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).
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).
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.

Regulation of inspiratory burst frequency and amplitude by I CAN /TRPM4 in preBötC glutamatergic neurons
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 . 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.

Model simulations-data comparisons
For comparisons of model predictions and results from the experimental photostimulation and pharmacological manipulations (Figure 8), we incorporated in the original model  a ChR2 current ( I ChR2 ) to represent I App for the excitatory neuron population and mimic effects of Figure 5. Power spectrum analyses of preBötzinger complex (preBötC) neuronal population activity in vitro before and after complete block of neuronal persistent sodium current (I NaP ). (A) Examples of preBötC integrated population activity patterns and associated power spectra before and after block of I NaP , which eliminates rhythmic preBötC population activity, leaving only a baseline low level of tonic activity or 'noise' with a flat power spectrum (lower traces). The power spectra of the rhythmic activity have clear peaks corresponding to the fundamental and higher harmonic frequencies (upper traces). (B) Representative experiment illustrating activity patterns and power spectra before and after block of I NaP at various levels of photostimulation (0 and 1 mW in control; 0, 1, 2, and 3 mW after I NaP blocked). With I NaP intact (control conditions, left traces), photostimulation increases the frequency of integrated preBötC inspiratory population activity accompanying an upward shift of the integrated baseline activity due to tonic activity, and there are clear peaks in the power spectra of the rhythmic activity corresponding to the fundamental and higher harmonic frequencies (lower panel). After block of I NaP with 20 nM tetrodotoxin citrate (TTX), there is a graded shift in the baseline level of tonic activity during photostimulation, but no rhythmic integrated population activity as indicated by the flat power spectra resembling that of the baseline noise activity in the absence of photostimulation. (C) Data from integrated preBötC population recordings showing the change in level of tonic activity relative to baseline during graded 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 photostimulation and the recovery period following stimulation under control conditions with I NaP intact and with I NaP blocked (20 nM TTX). Increasing the level of photostimulation (0.25, 1 mW in control; 1, 2, 3 mW with I NaP blocked) in both conditions increases the level of tonic activity. The activity returns to the baseline level (I NaP blocked) or is slightly depressed (control conditions) in the post-photostimulation recovery period. Bars indicate mean ± SEM from 5 slices. * indicates statistically significant difference from unity (p<0.05 by two-tailed t-test). (D) Quantification of the oscillatory component in the preBötC activity before and after complete block of I NaP at different levels of photostimulation as the ratio of the band power in the 0-3 Hz frequency range over the band power in the 3-6 Hz range. The ratios are high in control conditions, reflecting significantly higher spectral power content in the low frequency band (oscillations) compared to one in the high frequency band (noise) when there is rhythmic activity at baseline, with photostimulation, and during recovery. The ratios are not different and have a value of unity when I NaP is blocked, reflecting the flat power spectra under this condition. Bars indicate mean ± SEM from 5 slices.
The online version of this article includes the following source data for figure 5: Source data 1. Related to Figure 5C and D. , photostimulation did not induce rhythmic activity in this neuron, only membrane depolarization without rhythmic synaptic drive potentials or bursting at 1 mW laser application (B) and only tonic neuronal spiking at a higher laser power of 2 mW (C). The latter indicates that the neuron retained spiking capabilities at the low concentration of TTX employed, which did not interfere with action potential generation by transient Na + channels while I NaP was completely blocked (see Figure 4-figure supplement 1). (D) Voltageclamp recordings from td-tomato-labeled preBötC inspiratory neuron showing inward rhythmic synaptic drive currents synchronized with inspiratory hypoglossal activities. Under control conditions, optogenetic stimulation (1 mW) induced inward currents and rhythmic synaptic drive currents synchronized with the higher frequency inspiratory XII activities. (E) Rhythmic inspiratory synaptic drive currents did not occur with photostimulation (1 mW) after complete block of I NaP at 20 nM TTX, indicating loss of synaptic interactions and network rhythmic activity.
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 XII inspiratory activity (∫XII) during bath application of the specific pharmacological inhibitor of TRPM4 channels (9-phenanthrol, 50 µM), which gradually decreased inspiratory burst amplitude (~58% reduction at steady state ~60 min after drug application), but had little effect on the inspiratory frequency (bottom trace). (B) Bilateral photostimulation of the preBötC (control, before pharmacological block) increased inspiratory burst frequency in a laser-powerdependent manner (~284% increase at 1 mW, ~354% increase at 2 mW) and monotonically decreased inspiratory XII burst amplitude (~25% decrease at 1 mW, ~44% decrease at 2 mW). (C) Under TRPM4 block conditions (>60 min after 9-phenanthrol application), photostimulation significantly increased inspiratory frequency (~586% increase at 1 mW, ~724% increase at 2 mW), and decreased ∫XII inspiratory burst amplitude (~72% decrease at 1 mW, ~78% decrease at 2 mW). (D) Summary plots of monotonic relations between laser power and normalized inspiratory frequency (black lines, n=8 slices) and normalized amplitude (red, n=8). Data points plotted are mean ± SEM.
The online version of this article includes the following source data for figure 7: Source data 1. Related to Figure 7D. power curve is larger with riluzole than TTX application (as also shown in Phillips and Rubin, 2019) as seen in the experimental data.
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.

Model description and methods
As in the previous paper , 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 . In the updated model, the membrane potential Vm for each neuron is given by the following current balance equation: Cm dVm dt + I Na + I K + I Leak + I NaP + I CAN + I Ca + I Syn + I ChR2 = 0 , where Cm is the membrane capacitance, I Na , I K , I Leak , I NaP , I CAN , I Ca , I Syn , and I ChR2 are ionic currents through sodium, potassium, leak, persistent sodium, calcium-activated non-selective cation, voltagegated calcium, synaptic and ChR2 channels, respectively. For a full description of all currents other than I ChR2 , see . 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, δ is the ratio of the two open-state conductances, Vm is the membrane potential, and E ChR2 is the reversal potential for ChR2. All parameters are given in Table 1.
is an empirically derived voltagedependent 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 and burst frequency (top) and amplitude (bottom) of the preBötC network with model simulated blockade of I NaP by tetrodotoxin citrate (TTX) and riluzole (RZ) or block of I CAN . I NaP block by TTX and I CAN block by 9-phenanthrol (Phen.) were simulated by reducing the persistent sodium conductance (g NaP ) and calcium-activated non-selective cation conductance (g CAN ), respectively. In contrast, blockade of I NaP by RZ 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 RZ application, as described and justified previously (Phillips and Rubin, 2019). (C) Relationship between laser power and inspiratory burst frequency (top) and amplitude (bottom) of the integrated XII output with bath application of the pharmacological blockers of I NaP (TTX or RZ) or I CAN (9-phenanthrol); data are the same as shown in Figures 4C and 7D. Data points and error bars are mean ± SEM.
The online version of this article includes the following source data for figure 8: Source data 1. Related to Figure 8A-C.  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 ω 1 , ω 2 , γ 1 , γ 2 , γ 3 , γ 4 , and γ 5 are defined as follows: ω 1 = 0.8535 · F · P , ω 2 = 0.14 · F · P , γ 1 = 4.34 · 10 5 · e −0.0211539274·Vm , γ 2 = 0.008 + 0.004 · ln , T = 308.0 K , F = 96.485 kC/mol , , and γ 5 = 0.05 . 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 σret is the absorption cross-section, λ is the wavelength of max absorption, ω 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 τ ChR2 is the time constant of ChR2 activation. Model parameter values are given in Table 1.

Model tuning
Initial model predictions presented in Figure 1 were generated without modification of the initial model tuning presented in . 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.

Integration methods
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.

Animal procedures
All animal procedures were approved by the Animal Care and Use Committee of the National Institute of Neurological Disorders and Stroke (Animal Study Proposal number: 1154-21).

Cre-dependent triple transgenic mouse model for optogenetic manipulation of VgluT2-expressing glutamatergic neurons
We have previously established and histologically validated the Cre-dependent double transgenic mouse line (VgluT2-tdTomato) to study population-specific roles of VgluT2-expressing and red fluorescent protein-labeled glutamatergic preBötC neurons (Koizumi et al., 2016). This VgluT2-tdTomato mouse line was produced using a Slc17a6 tm2 (
We identified, under current-clamp recording, optically imaged (below) preBötC VgluT2-tdTomato expressing inspiratory neurons that have intrinsic voltage-dependent oscillatory bursting properties (endogenous burster neurons; Koizumi et al., 2013). As shown previously (Koizumi et al., 2013), these endogenous burster neurons exhibited 'ectopic' bursts that are not synchronized with inspiratory XII activity in addition to exhibiting inspiratory bursts synchronized with this inspiratory activity. Membrane depolarization of these neurons by applying steady current caused an increase of ectopic bursting frequency in a voltage-dependent manner as shown previously (Koizumi et al., 2013). In some experiments, we confirmed that these neurons exhibited intrinsic voltage-dependent rhythmic bursting behavior after blocking excitatory synaptic transmission by bath application of the non-NMDA glutamate receptor blocker, 6-cyano-7-nitroquinoxaline-2,3-dione disodium (CNQX) (20 µM; Sigma; also see Koizumi et al., 2013).
Whole-cell voltage-clamp recording from optically identified VgluT2-tdTomato expressing inspiratory neurons was used to obtain neuronal 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. Data acquisition and analyses were performed with PatchMaster and Igor Pro (Wavemetrics) software (Koizumi and Smith, 2008). Values of inspiratory neuronal membrane capacitance (C m ) and maximum I NaP conductance (g NaP ), were measured from voltage-clamp data as previously described (Koizumi and Smith, 2008).

Optical imaging of preBötC glutamatergic neurons
Optical two-photon live imaging of tdTomato-ChR2-EYFP-co-labeled neurons in the slice preparations was performed to verify fluorescent protein expression and target VgluT2-positive preBötC neurons for whole-cell recording. Images were obtained with a Leica multi-photon laser scanning upright microscope (TCS SP5 II MP with DM6000 CFS system, LAS AF software, 20× water-immersion objective, N.A. 1.0, Leica; and 560 nm beam splitter, emission filter 525/50, Semrock). A two-photon Ti:Sapphire pulsed laser (MaiTai, Spectra Physics, Mountain View, CA) was used at 800-880 nm with DeepSee predispersion compensation. The laser for two-photon fluorescence imaging was simultaneously used for transmission bright-field illumination to obtain a Dodt gradient contrast structural imaging to provide fluorescence and structural images matched to pixels. These images allowed us to accurately place a patch pipette on the tdTomato-labeled VgluT2-positive neurons to functionally identify preBötC inspiratory glutamatergic neurons that were active in phase with XII inspiratory activity.

Photostimulation of preBötC glutamatergic neuron population
Laser illumination for optogenetic experiments was performed with a blue laser (473 nm; OptoDuet Laser, IkeCool, Los Angeles, CA), and laser power (0.25-5.0 mW) was measured with an optical power and energy meter (PM100D, ThorLabs, Newton, NJ). Illumination epochs (sustained 20 Hz, 20 ms pulses, variable duration pulse trains) were controlled by a pulse stimulator (Master-8, A.M.P.I., Jerusalem, Israel). Optical fibers, from a bifurcated fiberoptic patch cable, each terminated by an optical cannula (100 µm diameter, ThorLabs), were positioned bilaterally on the surface of the preBötC in the in vitro rhythmically active slice preparations.

Signal analysis of respiratory parameters and statistics
All digitized electrophysiological signals were analyzed off-line to extract respiratory parameters (inspiratory burst frequency and amplitude) from the smoothed integrated XII nerve or preBötC neuronal population inspiratory activities (200ms window moving average) performed with Chart software (AD Instruments), Spike2 software (Cambridge Electronics Design), and Igor Pro (WaveMetrics) software. Following automated peak detection, burst period was computed to obtain the inspiratory burst frequency. Inspiratory burst amplitude was calculated as the difference between the signal peak value and the local baseline value. To obtain steady-state perturbations, we analyzed the respiratory parameters during laser illumination excluding the first and last 5-10 s of laser illumination. Power spectrum analyses were performed using MATLAB software ver. R2020b (MathWorks, Natick, MA, USA). Power spectra were calculated as a squared magnitude of the Fast Fourier Transform of the rectified preBötC activity for 30 s segments of the recordings during baseline activity, photostimulation, and recovery periods. The oscillatory component in preBötC activity was quantified as a ratio of the band power of the signal in the 0-3 Hz frequency range to the band power of the signal in the 3-6 Hz frequency range. For statistical analyses, respiratory parameters during laser illumination were normalized to the control value calculated as an average from ~30 inspiratory bursts before laser illumination. The normality of data was assessed both visually (quantile vs quantile plots) and through the Shapiro-Wilk normality test. Statistical significance (p<0.05) was determined with non-parametric Wilcoxon matched-pairs signed rank test or Kolmogorov-Smirnov test when comparing two groups, and two-way ANOVA test for comparing multiple groups in conjunction with post hoc Tukey's HSD test for pairwise comparisons (Prism, GraphPad software LLC), and summary data are presented as the mean ± SEM.

Discussion
In these new combined experimental and modeling studies, we have advanced 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. These studies were designed to further experimentally test predictions of our recent computational model  of preBötC excitatory circuits that incorporate rhythm-and amplitude-generating biophysical mechanisms, relying, respectively, on a neuronal subthreshold-activating, slowly inactivating I NaP , and an I CAN , mediated by TRPM4 channels coupled to intracellular calcium dynamics. Our model explains how the functions of generating the rhythm and amplitude of inspiratory oscillations in preBötC excitatory circuits involve distinct biophysical mechanisms. In essence, the model advances the concepts that (1) a subset of excitatory circuit neurons whose rhythmic bursting in vitro is critically dependent on I NaP forms an excitatory neuronal kernel for inspiratory rhythm generation, and (2) excitatory synaptic drive from the rhythmogenic kernel population is critically amplified by I CAN activation in recruited and interconnected preBötC excitatory neurons in the network to generate the amplitude of population activity.
The model predicts that I NaP is critical for rhythm generation and confers voltage-dependent rhythmic burst frequency control in vitro, which provides a mechanism for frequency tuning over a wide dynamic range defined by the frequency tuning curve for the rhythmogenic population (i.e. the relationship between applied current/baseline membrane potential across the network and network bursting frequency). Accordingly, reducing neuronal g NaP decreases the population-level bursting frequency with a weaker reduction in amplitude, reduces the frequency dynamic range, and, at sufficiently low levels of g NaP , terminates rhythm generation under baseline and other conditions of network excitability in vitro, as also shown in Phillips and Rubin, 2019. In contrast, I CAN in the model is essential for generating the amplitude of rhythmic output but not for rhythm generation. As a result, reducing g CAN strongly decreases population activity amplitude and has little effect on population bursting frequency under baseline conditions in vitro, but strongly augments bursting frequency, due to a shift in the frequency tuning curve that extends the frequency range, at higher levels of population excitation. These predicted opposing effects of I NaP and I CAN attenuation on the relationship between network excitability and preBötC population rhythmic activity provided a clear basis for model testing, and our experimental results showed an overall strong agreement with the model predictions.
To experimentally test these predictions in the rhythmically active in vitro preBötC network in slices from transgenic mice, we used a combination of electrophysiological analyses, pharmacological perturbations of I NaP or I CAN , and selective optogenetic manipulations of the preBötC excitatory (glutamatergic) population involved. Our application of optogenetic photostimulation of preBötC glutamatergic neurons to examine the population activity frequency tuning curves under these different conditions of I NaP or I CAN attenuation provided an effective way to probe for underlying mechanisms. We performed new simulations to mimic the optogenetic manipulations under conditions of pharmacological perturbations of I NaP or I CAN . I NaP and I CAN block in the experimental data yield shifts in the preBötC network burst frequency/amplitude vs laser power curves (Figure 8) in the same direction, and of the same relative magnitudes, as predicted by our simulations (Figures 1 and 8). This correspondence strongly supports the main hypothesis from the model that I NaP and I CAN mechanistically underlie preBötC rhythm and inspiratory burst pattern generation, respectively. The similarities between the model predictions and experimental results support the basic concepts about inspiratory rhythm and burst pattern generation from the model.

Extensions and limitations of the model
The assumptions and neurobiological simplifications of the original model were previously elaborated in  We have extended the model to allow comparisons of our new experimental data and model simulation results. The updates included incorporating a I ChR2 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 , although the details of the voltage-dependent kinetic properties are more complex than represented by the firstorder 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  Further insights into mechanisms of rhythm and burst amplitude generation in preBötC excitatory circuits in vitro 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;. 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  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 2+ 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.
Previous pharmacological studies and proposed roles of I NaP in preBötC inspiratory network rhythm generation 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 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 and Rubin, 2019) and experimental studies (Johnson et al., 1994). Moreover,  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 2+ imaging in glutamatergic neurons in vitro .

Summary
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.

Additional files
Supplementary files • Transparent reporting form • Source code 1. Model source code.

Data availability
All data generated or analyzed during this study are included in the manuscript and supporting source data files. All model parameters and equations are included in the manuscript and source code is provided.