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. 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; (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 neonatal mouse slices in vitro. 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.


1
Model concepts and predictions 2 The following features of the model outputs from new simulations (Figure 1 (1) The inspiratory rhythm in vitro is generated by a functionally distinct subpopulation of 7 preBötC excitatory neurons with subthreshold-activating INaP that confers voltage-8 dependent rhythmic burst frequency control, providing a mechanism for frequency tuning 9 by the regulation of baseline membrane potentials. Progressively 10 depolarizing/hyperpolarizing neurons across the population by varying an applied current 11 (IApp) increases/decreases population burst frequency over a wide dynamic range defined 10nM TTX, n=6/7 at 5 µM riluzole, n=8/8 at 10 µM riluzole). In the cases where rhythmic 1 bursting was terminated, after 30 min. we again bilaterally photostimulated the preBötC, 2 which restored rhythmic bursting with voltage-dependent frequency control (0.25 -2.0 mW) 3 in all cases (n=14/14 with TTX, n=14/14 with riluzole). The summary data (n=6 slice 4 preparations for TTX 5nM, n=8 for TTX 10 nM, n=6 for riluzole 5 µM, n=8 for riluzole 10 µM, 5 mean ± SEM) in Figure 4C (left) shows the frequency tuning curves with significant increases 6 in burst frequency as a function of laser power (one sample t-test, p<0.01 in all cases except 7 p<0.05 at 0.25 mW with 10 nM TTX and at 0.25 mW with 10 µM riluzole). The tuning curves 8 are significantly shifted downward relative to control conditions (TTX 5 nM vs control, TTX 9 10 nM vs control, riluzole 5 µM vs control, riluzole 10 µM vs control; in all cases, p<0.01 by 10 Wilcoxon matched pair signed rank test), and more downwardly-shifted with higher drug 11 concentrations (Wilcoxon signed rank test, TTX 5 nM vs TTX 10 nM, p<0.01; riluzole 5 µM vs 12 riluzole 10 µM, p = 0.018). The frequency-laser power relations are not significantly different 13 for TTX and riluzole (Wilcoxon signed rank test, p = 0.101 for TTX 5 nM vs riluzole 5 µM; p = 14 0.431 for TTX 10 nM vs riluzole 10 µM). 15 Summary data (n = 6 slice preparations at TTX 5nM, n=8 at TTX 10 nM, n=6 at riluzole 5 16 µM, n=8 at riluzole 10 µM, mean ± SEM) of the relations between normalized inspiratory 17 burst amplitude and laser power presented in Figure 4C  significantly downward compared to control conditions (TTX 5 nM vs control, p = 0.290; TTX 23 10 nM vs control, p= 0.036; riluzole 5 µM vs control, p<0.01; riluzole 10 µM vs control, p<0.01 1 by Wilcoxon signed rank tests), but there were no significant differences (Wilcoxon signed 2 signed rank test, p<0.01). Note that under baseline conditions of excitability, TRPM4 block 1 does not alter burst frequency, and the upward shift in the frequency tuning curve is 2 opposite to the shift seen with TTX or RZ block of INaP. These findings are consistent with 3 model predictions presented in Figure 1.  TTX and riluzole on INaP, there is a larger reduction of amplitude with INaP block by riluzole than by TTX (as also shown in Phillips and Rubin,  1 PLoS CB, 2019), as seen in the experimental data. 2 Reducing gCAN more strongly decreases population activity amplitude in the model 3 simulations and, in contrast to reducing gNaP, has little effect on population bursting 4 frequency under baseline conditions, but strongly augments bursting frequency at higher 5 levels of excitability. This results in an upward shift in the frequency tuning curve that 6 extends the frequency range at the higher levels of population depolarization, consistent 7 with the experimental results ( Figure 7). Thus, the model predictions are qualitatively 8 consistent with the experimental data for major features of excitatory network behavior. 9

10
The rates of change between states are given by the following differential equations: 11 12 4 = 4 $ + 7 4 − 4 4 13 irradiance-dependent activation function for ChR2, respectively. These quantities take the form 15 where A-B is the absorption cross-section, is the wavelength of max absorption, CDEE a scaling 1 factor for loss of photons due to scattering/absorption, ℎ is Planks constant, is the speed of light, 2 and %&'$ is the time constant of ChR2 activation. Model parameter values are given in Table 1. 3 4

Model Tuning 5
In order to match experimental data, parameters were slightly modified from the original model 6 tuning in Phillips et al. (2019). The updated list of parameters is given in Table 1, where parameter 7 adjustments are indicated with red font.

Integration methods 2
All simulations were performed locally on an 8-core Linux-based operating system. Simulation 3 software was custom written in C++. Numerical integration was performed using the exponential We recorded motor outputs in vitro from XII nerve rootlets with fire-polished glass suction 5 electrodes (50 -100 µm inner diameter) to monitor inspiratory XII motoneuron population 6 activity. In addition, to monitor preBötC inspiratory population activity directly, we 7 performed macro-patch recordings by applying a fire-polished glass pipette (150 -300 µm 8 inner diameter) directory on the surface of the preBötC region. In some experiments, 9 extracellular neuronal activity was also recorded from preBötC inspiratory neurons with a 10 fine-tipped glass pipette (3 -5 MΩ resistance) filled with 0.5 M sodium acetate (Sigma) 11 positioned by a computer-controlled 3-dimensional micromanipulator (MC2000, gluconate, 5.0 Na-gluconate, 3.0 NaCl, 10.0 HEPES, 4.0 Mg-ATP, 0.3 Na-GTP, and 4.0 sodium 1 phosphocreatine, pH 7.3 adjusted with KOH. In all cases, measured potentials were corrected 2 for the liquid junction potential (-10 mV). Series resistance was compensated on-line by 3 ~80%, and the compensation was periodically readjusted. 4 Whole-cell voltage-clamp recording from optically identified VgluT2-tdTomato expressing 5 inspiratory neurons was used to obtain subthreshold current-voltage (I-V) relationships by 6 applying slow voltage ramps (30 mV/s; -100 to +10 mV) and we computed TTX-and riluzole- compensation. The laser for two-photon fluorescence imaging was simultaneously used for 20 transmission bright-field illumination to obtain a Dodt gradient contrast structural imaging 21 to provide fluorescence and structural images matched to pixels. These images allowed us to 22 accurately place a patch pipette on the tdTomato-labeled VgluT2-positive neurons to functionally identify preBötC inspiratory glutamatergic neurons, which were active in phase 1 with XII inspiratory network activity. Electronics Design), and Igor Pro (WaveMetrics) software. Following automated peak 19 detection, burst period was computed to obtain the inspiratory frequency. Inspiratory burst 20 amplitude was calculated as the difference between the signal peak value and the local 21 baseline value. To obtain steady-state perturbations, we analyzed the respiratory 22 parameters during laser illumination excluding the first and last 5-10 sec of laser illumination. Power spectrum analyses were performed using Matlab software. Power 1 spectra were calculated as a squared magnitude of the Fast Fourier Transform of the rectified 2 preBötC activity for 30 sec segments of the recordings during baseline activity, 3 photostimulation and recovery period. The oscillatory component in preBötC activity was 4 quantified as a ratio of the band power of the signal in the 0 -3 Hz frequency range to the 5 band power of the signal in the 3 -6 Hz frequency range. For statistical analyses, respiratory 6 parameters during laser illumination were normalized to the control value calculated as an 7 average from ~30 inspiratory bursts before laser illumination. The normality of data was 8 assessed both visually (Quantile vs Quantile Plots) and through the Shapiro-Wilk normality 9 test. Statistical significance (p < 0.01 or p < 0.05) was determined with Student's t-test or 10 Wilcoxon signed rank test (Prism, GraphPad software LLC), and summary data are presented 11 as the mean ± SEM. The model predicts that INaP is critical for rhythm generation and confers voltage-dependent 17 rhythmic burst frequency control, which provides a mechanism for frequency tuning over a wide 18 dynamic range defined by the frequency tuning curve for the rhythmogenic population (i.e., the 19 relationship between applied current/baseline membrane potential across the network and network 20 bursting frequency). Accordingly, reducing neuronal persistent sodium conductance (gNaP) 21 decreases the population-level bursting frequency with a weak reduction in amplitude, reduces the 22 frequency dynamic range, and, at sufficiently low levels of gNaP, terminates rhythm generation under baseline and other conditions of network excitability in vitro, as also shown in Phillips and 1 Rubin, 2019. In contrast, ICAN in the model is essential for generating the amplitude of rhythmic 2 output but not for rhythm generation. As a result, reducing gCAN strongly decreases population 3 activity amplitude and has little effect on population bursting frequency under baseline conditions 4 in vitro, but strongly augments bursting frequency, due to a shift in the frequency tuning curve that 5 extends the frequency range, at higher levels of population excitation. These predicted opposing 6 effects of INaP and ICAN attenuation on the relationship between network excitability and preBötC 7 population rhythmic activity provided a clear basis for model testing, and our experimental results 8 showed an overall strong agreement with the model predictions. 9 To experimentally test these predictions in the rhythmically active in vitro preBötC network in 10 slices from transgenic mice, we used a combination of electrophysiological analyses,  (Williams et al., 2013). We found that this 7 biophysical model could be parameterized to closely match our experimental data on relations 8 between laser power and membrane depolarization of identified preBötC inspiratory glutamatergic 9 neurons. Furthermore, we found that this channel model parameterized from our cellular-level data 10 yielded frequency tuning curves for the bilateral excitatory population under baseline conditions 11 and with pharmacological perturbations that were in close agreement with the experimentally 12 measured tuning curves. We note that in the model, the assumption is that all excitatory network burst-terminating mechanisms incorporated in these models are apparently not sufficiently 5 expressed in the in vitro preBötC network. 6 We note that our new measurements indicate that while I !UV /TRPM4 activation is not involved 7 in rhythm generation under baseline conditions of network excitability, since inhibiting this current 8 does not significantly affect the rhythm, we also show that reducing this current does shift the 9 frequency tuning curve at higher levels of network excitation. This occurs because activation of 10 this current augments excitatory synaptic interactions, which can affect network bursting 11 frequency. As discussed above, these observed effects of reducing I !UV /TRPM4 on bursting 12 frequency at elevated levels of network excitation, which we were able to reveal with our 13 photostimulation paradigm, is a basic feature and entirely consistent with predictions of our model. 14 Our new experimental results confirm our previous experimental and modeling results 15 indicating that endogenous activation of I !UV /TRPM4 is critically involved in generating the 16 amplitude of population activity. The correspondence between the experimental data and model 17 predictions supports the concept in the model that activation of ICAN is largely due to synaptically-18 activated sources of neuronal Ca 2+ flux such that ICAN contributes to the excitatory inspiratory drive 19 potential and regulates inspiratory burst amplitude by augmenting the excitatory synaptic current. 20 Our data also show that INaP is involved, to a small degree, in generating the amplitude of rhythmic 21 population activity at a given level of excitatory network excitation. This is due to the subthreshold 22 activation of INaP and its voltage-dependent amplification of synaptic drive; indeed, application of riluzole, which is thought to impact synaptic transmission, affected population amplitude much 1 more than application of TTX (Figure 4, Figure 8). This contribution to amplitude is smaller than 2 that from I !UV /TRPM4 activation, however. In general, our results support the concept from our 3 original model that I !UV activation in a subpopulation of preBötC excitatory neurons is critically 4 involved in amplifying synaptic drive from a subset of neurons whose rhythmic bursting is 5 critically dependent on )*/ . This later subpopulation forms the kernel for rhythm generation in In addition, thicker slices than used in our studies contain not only the preBötC, but also other 11 network components that may introduce additional mechanisms, which is why we use thin slice 12 preparations to allow intrinsic properties of the isolated but functional preBötC circuits to be more 13 Regardless of this greater complexity, our present study has confirmed our previous results that 1 INaP, is essential for inspiratory rhythm generation in the isolated preBötC excitatory circuits in 2 vitro, representing a major rhythmogenic mechanism in this reduced state of the respiratory 3 network. An important prediction of our model (see Phillips et al. (2019) Fig. 9) is that the 4 subpopulation of excitatory neurons with INaP-dependent rhythmogenic properties forming the 5 rhythmogenic kernel may be relatively small compared to the neuronal subpopulation(s) with the 6 ICAN-dependent mechanism that critically amplifies the rhythmic drive from the kernel to generate 7 and control the amplitude of excitatory population activity. This prediction also requires direct Under these conditions of partial block of INaP, bilateral preBötC photostimulation (30 min. after 7 the rhythm stopped, blue shaded epoch) could reinitiate the rhythm, which was also laser power 8 dependent (~202% increase at 1 mW compared to the control, and ~306% increase at 2 mW in the 9 example shown). (C) Summary plots (TTX 5nM; n=6, TTX 10 nM; n=8, riluzole (RZ) 5µM; n=6, 10 RZ 10 µM; n=8, mean ± SEM) of the relations between laser power vs normalized inspiratory 11 burst frequency indicates laser power dependent, significant increases of frequency in all cases. 12 The curves for higher concentration of TTX and riluzole are downward-shifted compared to those 13 for the lower concentration as well as those under control conditions (before drug applications).
Right panel shows summary plots (TTX 5nM; n=6, TTX 10 nM; n=8, RZ 5µM; n=6, RZ 10 µM; 1 n=8, mean ± SEM) of the relations between laser power vs normalized XII burst amplitude 2 indicating laser power dependent, significant decreases of burst amplitude in all cases. The curves 3 under the INaP partial block (except for at TTX 5 nM) are downward-shifted compared to those 4 under the control conditions. The decrease in burst amplitudes under RZ is more significant than 5 those under TTX, although there are no significant differences between different concentrations 6 of TTX (blue lines) or RZ (red lines).