Theory of hierarchically organized neuronal oscillator dynamics that mediate rodent rhythmic whisking

Abstract

Rodents explore their environment through coordinated orofacial motor actions, including whisking. Whisking can free-run via an oscillator of inhibitory neurons in the medulla and can be paced by breathing. Yet, the mechanics of the whisking oscillator and its interaction with breathing remain to be understood. We formulate and solve a hierarchical model of the whisking circuit. The first whisk within a breathing cycle is generated by inhalation, which resets a vibrissa oscillator circuit, while subsequent whisks are derived from the oscillator circuit. Our model posits, consistent with experiment, that there are two subpopulations of oscillator neurons. Stronger connections between the subpopulations support rhythmicity, while connections within each subpopulation induce variable spike timing that enhances the dynamic range of rhythm generation. Calculated cycle-to-cycle changes in whisking are consistent with experiment. Our model provides a computational framework to support longstanding observations of concurrent autonomous and driven rhythmic motor actions that comprise behaviors.

Keywords: brainstem; breathing; medulla; network; oscillations; rate modeling; synchrony; vibrissa.

Figure 1.. Analysis of existing and new data related to breathing and whisking.

(A) Typical experimental set-up for recording from head-restrained rodents. (B,C) Time series of breathing (red) and whisking (blue) for a head-restrained rat [Moore et al. (2013)]. (D) Annotated whisking and breathing. The values ${\theta }_1^{\text{amp}},{\theta }_2^{\text{amp}}$, and ${\theta }_3^{\text{amp}}$ denote the amplitudes of the first, second and third whisks within the same breathing cycle. (E) Time series of breathing (red) and whisking (blue) for a free moving rat. (F) Semi-logarithmic plot of the number of breaths for which there are $N_{\mathrm{W}}$ whisking cycles. The dotted lines are exponential fits with slopes of 0.95 (restrained) and 0.92 (free). (G) Amplitude of the first whisk ${\theta }_1^{\text{amp}}$ within a breathing cycle versus $N_{\mathrm{W}}$. Bars denote the standard error. (H) Probability density function of $\mathrm{l}\mathrm{o}\mathrm{g}\left({\theta }_{I+1}^{\mathrm{a}\mathrm{m}\mathrm{p}}/{\theta }_i^{\text{amp}}\right)$ versus ${\theta }_i^{\text{amp}}$ for $i=1,2,3$. The normalization of each distribution is ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }} d\left[\mathrm{l}\mathrm{o}\mathrm{g}\left({\theta }_{I+1}^{\mathrm{a}\mathrm{m}\mathrm{p}}/{\theta }_i^{\mathrm{a}\mathrm{m}\mathrm{p}}\right)\right]\mathrm{l}\mathrm{o}\mathrm{g}\left({\theta }_{I+1}^{\mathrm{a}\mathrm{m}\mathrm{p}}/{\theta }_i^{\mathrm{a}\mathrm{m}\mathrm{p}}\right)=1$. (I) Modulation depth of the spike rates of vIRt neurons. Black circles are from intracellular recordings with lightly anesthetized rats [Deschênes et al. (2016b)]. Green circles represent new data from awake mice with extracellular electrodes and optical tagging. The depth is calculated from the spike rates as $<($ minimum – maximum $)/$ average $>$ where the averaging is on a per cycle basis. Data plotted as a function of the phase in the whisk cycle at which the spike rate is maximal. (J) New analysis to determine the coefficient of variation ${\mathrm{C}\mathrm{V}}_2$ as a function of the phase in the whisk cycle. Same data and notations as in panel I.

Figure 2.. Architecture of the brainstem circuit model for whisking.

(A) The neuronal-level circuit for conductance-based modeling. The triangles are neurons. The blue and red colors denote inhibitory and excitatory connections, respectively. The currents $I_{\mathrm{e}\mathrm{x}\mathrm{t}}^{\mathrm{r}}$ and $I_{\mathrm{e}\mathrm{x}\mathrm{t}}^{\mathrm{F}}$ represent constant external depolarizing input to the vIRt and vFNM neuronal populations, respectively. Conductances between neuronal pairs are denoted by $g_{\text{intra}}/K$ for pairs of neurons that belong to the same subpopulation, by $g_{\text{inter}}/K$ for neuronal pairs from two different subpopulations, and by $g_{\mathrm{F}\mathrm{r}}/K$ for vIRt^ret-to-vFNM connections. The amplitude of the square-wave pBötC-to-vIRt^ret input is denoted by $g^{\mathrm{r}\mathrm{B}}$. (B) Schematic of the different currents in the cellular modelfor vIRt cells; the same model applies to vFM neurons with possibly different conductances. The currents are a leak current, ${\mathrm{I}}_{\mathrm{L}}$, the transient sodium current, $I_{\mathrm{N}\mathrm{a}}$, the delayed rectifier potassium current, $I_{\mathrm{K}\mathrm{d}\mathrm{r}}$, the persistent sodium current $I_{\mathrm{N}\mathrm{a}\mathrm{P}}$, the mixed cation h-current, $I_{\mathrm{h}}$, an M-type ${\mathrm{K}}^{+}$ current, $I_{\text{adapt}}$, external excitatory currents from other brain areas, $I_{\text{ext}}$ and synaptic currents $I_{\mathrm{s}\mathrm{y}\mathrm{n}}$ that comprise $I_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}},I_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$ and $I^{\mathrm{r}\mathrm{B}}$ or $I^{\mathrm{F}\mathrm{r}}$. Details in Star Methods. (C) Spiking rate versus $I_{\text{ext}}^{\mathrm{r}}$ curves of the single vIRt neuron model for values of $g_{\text{adapt}}$ that range from 0 to 7 mS/cm² (solid lines). The insert is for $g_{\text{adapt}}=6\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2$ with the spike rate versus $I_{\mathrm{e}\mathrm{x}\mathrm{t}}$ curve as solid and a linear approximation as dotted, as used in the rate model. (D) The force developed by an intrinsic muscle that moves vibrissae as a function of the average firing rate of vFMN neurons, $M^{\mathrm{F}}$. The solid black line denotes simulation results, and the dotted red line denotes a fit. Details in Star Methods. (E) Schematic of the rate model in which each population of neurons is described by an activation variable that controls its synaptic outputs, and another activation variable that controls the adaptation current. The interactions $J_{\text{intra}}$ and $J_{\text{inter}}$ replace $g_{\text{intra}}$ and $g_{\text{inter}}$, respectively.

Figure 3.. Dynamics of a conductance-based and rate-based feed-forward circuit for whisking.

(A) Schematic of the conductance-based circuit with input from the pBötC to the vIRt^ret and from the vIRt^ret to the vFMN to drive whisking. The vIRt^pro subpopulation plays no role. We further set $g_{\text{intra}}=0$. We chose $g^{\mathrm{r}\mathrm{B}}=0.5\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2,g_{\mathrm{F}\mathrm{r}}=3\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2,I_{\text{ext}}^{\mathrm{r}}=20\mu \mathrm{A}/{\mathrm{c}\mathrm{m}}^2,I_{\text{ext}}^{\mathrm{F}}=3.1\mu \mathrm{A}/{\mathrm{c}\mathrm{m}}^2,K=25,g_{\text{adapt}}^{\mathrm{r}}=7\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2,g_{\text{adapt}}^{\mathrm{F}}=0.3\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2,N=100,T_{\mathrm{p}\mathrm{B}\mathrm{ö}\mathrm{t}\mathrm{C}}=200\mathrm{m}\mathrm{s},T_{\text{rand}}=10\mathrm{m}\mathrm{s},\mathrm{\Delta }{T}_{\mathrm{v}\mathrm{I}\mathrm{R}\mathrm{t}}=70\mathrm{m}\mathrm{s}$. Details of all calculations are in Star Methods. (B) Time courses of spiking for an example neuron from the vIRt^ret and vFMN subpopulations. (C) Rastergrams across twenty neurons in the simulated subpopulations. (D) The calculated vibrissa angle $\theta \left(t\right)$, calculated with the model of the motor plant and common parameters. (E) The reduced rate-based circuit with input from the pBötC to the vIRt^ret and from the vIRt^ret to the vFMN to drive whisking. (F-G) Circuit properties as the synaptic conductance from the pBötC is increased. Parameters are as in panel A, except that $g^{\mathrm{r}\mathrm{B}}$ is varied. Properties are computed using three modeling strategies. First, by numerical simulations of the conductance-based model (solid black line and circles). Simulations were carried out over five realizations for each parameter set, and the error bars denote standard deviation. Second, by numerical simulations of the rate model equations (thick cyan line). Third, from analytical solution of the rate model (thin orange line). The average number of spikes produced by the vIRt^ret neurons per whisk is shown in panel F. The average number of spikes produced by vFMN neurons per whisk is shown in panel G. The average amplitude of each whisk is shown in panel H.

Figure 4.. Dynamics of conductance-based circuits without pBötC input to the vIRt.

(A) Schematic of the circuit with connections only between vIIRt^ret and vIRt^pro subpopulations, for which we set $g_{\text{intra}}=0$. We chose $g_{\text{inter}}=6\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2,g_{\text{adapt}}=7\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2,I_{\mathrm{e}\mathrm{x}\mathrm{t}}=20\mu \mathrm{A}/{\mathrm{c}\mathrm{m}}^2,K=25$, and $N=100$. (B) Time courses of spiking for an example neuron from the vIIRt^ret, vIIRt^pro, and vFMN subpopulations. (C) Rastergrams across twenty neurons in the simulated subpopulations. (D) The calculated vibrissa angle $\theta \left(t\right)$ with the model of the motor plant. (E) Schematic of the circuit with connections both between vIRt^ret and vIRt^pro subpopulations and within each subpopulation, for which we set $g_{\text{intra}}=12\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2$ and $g_{\text{inter}}=20\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2$; other parameters as in panel A (F) Time courses of spiking for an example neurons for the vIRt^ret, vIRt^pro, and vFMN subpopulations. Properties are computed using three modeling strategies as in Figure 3F–G. (G) Rastergrams across twenty neurons in the simulated subpopulations. (H) The calculated vibrissa angle $\theta \left(t\right)$.

Figure 5.. Dynamical properties of circuits without pBötC input to the vIRt.

(A) Schematic of the circuit. The dynamics are calculated for $g_{\text{intra}}=12\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2,g_{\text{adapt}}^{\mathrm{r}}=7\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2,g_{\text{adapt}}^{\mathrm{F}}=0.3\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2$ and $I_{\text{ext}}^{\mathrm{F}}=3.1\mu \mathrm{A}/{\mathrm{c}\mathrm{m}}^2$. The excitatory input fixed at $I_{\text{ext}}^{\mathrm{r}}=20\mu \mathrm{A}/{\mathrm{c}\mathrm{m}}^2$ and $g_{\text{inter}}$ is varied in terms of $\mathrm{\Delta }{g}_{\mathrm{s}\mathrm{y}\mathrm{n}}=g_{\text{inter}}-g_{\text{intra}}$ for panels B, D, F, H, and J). $\mathrm{\Delta }{g}_{\mathrm{s}\mathrm{y}\mathrm{n}}$ is fixed at $\mathrm{\Delta }{g}_{\mathrm{s}\mathrm{y}\mathrm{n}}=8\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2$ (black arrow in in panel L) with $I_{\text{ext}}$ varying for panels C, E, G, I, and K. Properties are computed using three modeling strategies as in Figure 3F–G, and we use the same notation. (B,C) The average spike rate $⟨M{⟩}_i$. The values for the vIRt^ret and vIRt^pro are equal in the uniform state, as the neuronal subpopulations are constantly spiking, and in the symmetric oscillatory state, since the neuronal subpopulations are alternately active. As of $\mathrm{\Delta }{g}_{\mathrm{s}\mathrm{y}\mathrm{n}}$ increases, one subpopulation becomes more active than the other and the less active subpopulation becomes silent at large values of $\mathrm{\Delta }{g}_{\mathrm{s}\mathrm{y}\mathrm{n}}$. The rate model exhibits a transition from a symmetric oscillatory state to a bistable state at a value of $\mathrm{\Delta }{g}_{\mathrm{s}\mathrm{y}\mathrm{n}}=g_{\mathrm{d}\mathrm{e}\mathrm{t}}$ while the actual transition in the conductance-based model occurs when of $\mathrm{\Delta }{g}_{\mathrm{s}\mathrm{y}\mathrm{n}}$ has further increased. The values of $g_{\mathrm{t}\mathrm{r}}$ and $g_{\mathrm{d}\mathrm{e}\mathrm{t}}$ are defined in Star Methods. (D,E) The whisking frequency $1/T_{\text{vIRt}}$. (F,G) The average whisking amplitude $<{\theta }^{\text{amp}}{>}_t$. Analytical results are computed in panel G for the uniform and the bistable state. (H,I) The whisking set-point ${<{\theta }^{\text{set}}>}_t$. Analytical results are computed in panel I for the uniform and the bistable state. (J,K) The coefficient of variation ${⟨{\mathrm{C}\mathrm{V}}_2⟩}_i$ calculated solely from the conductance-based equations. (L) Phase diagram showing the three dynamical regimes, uniform (left), oscillatory (middle) and bistable (right) computed using the conductance-based model. Values of the coefficient of variation ${⟨{\mathrm{C}\mathrm{V}}_2⟩}_i$ calculated for several values of $\mathrm{\Delta }{g}_{\text{syn}}$ and $g_{\text{intra}}$ are written. (M) The coefficient of variation ${⟨{\mathrm{C}\mathrm{V}}_2⟩}_i$ calculated as a function of K. The grey ribbon denotes typical experimentally-measured values for ${\mathrm{C}\mathrm{V}}_2$ (Figure 1J)

Figure 6.. Dynamics of a conductance-based circuits with pBötC input to the vIRt.

See Figure 2A for the schematic. We used the same parameters as the simulation without pBötC input (Figure 4E) plus $g^{\mathrm{r}\mathrm{B}}=0.5\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2,T_{\mathrm{p}\mathrm{B}\mathrm{ö}\mathrm{t}\mathrm{C}}=700\mathrm{m}\mathrm{s},T_{\text{rand}}=150\mathrm{m}\mathrm{s}$, and $\mathrm{\Delta }{t}_{\mathrm{p}\mathrm{B}\mathrm{ö}\mathrm{t}\mathrm{C}}=70\mathrm{m}\mathrm{s}$. (A) Time courses of spiking for an example neurons for the vIRt^ret, vIRt^pro, and vFMN subpopulations. (B) Rastergrams across twenty neurons in the simulated subpopulations. (C) The calculated vibrissa angle $\theta \left(t\right)$. (D) The calculated whisking amplitude as a function of $g_{\text{inter}}$. The first whisk is driven by the pBötC input and subsequent intervening whisks, with $i=2\text{-}4$, are driven by internal vIRt dynamics. Error bars in D-F denote SD. (E). The calculated whisking amplitudes as a function of $g^{\mathrm{r}\mathrm{B}}$. (F). The slowing down of whisking for the intermediate whisks is shown by plotting ${<t_{w,i+1}-t_{w,i}>}_t$ as a function of $g_{\text{inter}}$.

Figure 7.. Effects of pBötC activity on the timing and amplitude of the subsequent whisk.

(A). Cartoon to define the symbols. (B). Stimulations of the conductance model to show the vibrissa angle $\theta \left(t\right)$ for two aspects of the timing of breathing relative to whisking (Figure 4). The common vertical dotted line, labeled $t_{\mathrm{w},1}$, denotes the peak of a whisk just prior to input from the pBötC. The top trace is an example of input from the pBötC just after the time of the peak and the bottom trace is an example of relatively late input. Model parameters: $g^{\mathrm{r}\mathrm{B}}=0.5\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2,T_{\mathrm{p}\mathrm{B}\mathrm{ö}\mathrm{t}\mathrm{C}}=700\mathrm{m}\mathrm{s},T_{\text{rand}}=150\mathrm{m}\mathrm{s}$ (defined in Star Methods), $\mathrm{\Delta }{t}_{\mathrm{p}\mathrm{B}\mathrm{ö}\mathrm{t}\mathrm{C}}=70\mathrm{m}\mathrm{s}$. (C). Experimental results from head-restrained rats [Moore et al. (2013)] of the total time between protraction, $\mathrm{\Delta }{t}_{\mathrm{w},21}$, as a function of the time of inhalation after the last protraction, $\mathrm{\Delta }{t}_{\mathrm{B}\mathrm{w},1}$. The gray line is the linear fit over the full range. The mean is over all trials. (D). The calculated total time $\mathrm{\Delta }{t}_{\mathrm{w},21}$ as a function of $\mathrm{\Delta }{t}_{\mathrm{B}\mathrm{w},1}$ for two values of input conductance from the pBötC, i.e., $g^{\mathrm{r}\mathrm{B}}=0.5\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2$ (top) and $0.05\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2$ (bottom). The gray line is the linear fit over the range $40\mathrm{m}\mathrm{s}<\mathrm{\Delta }{t}_{\mathrm{B}\mathrm{w},1}<140\mathrm{m}\mathrm{s}$. For $g^{\mathrm{r}\mathrm{B}}=0.05\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2$ and small values of $\mathrm{\Delta }{t}_{\mathrm{B}\mathrm{w},1}$, the deduced values of $\mathrm{\Delta }{t}_{\mathrm{w},21}$ may be large since the first increase in the vibrissa angle caused by input from the pBötC may be too small to be detected as a separate new whisk. (E). The slope $s_{\mathrm{w},21}$ (top) and the intercept $\mathrm{\Delta }{t}_{\mathrm{w},21}\left(0\right)$ (bottom) as a function of $g^{\mathrm{r}\mathrm{B}}$. Error bars denote SD. The dotted lines denote the experimental values from panel D. (F). The whisking amplitude ${<{\theta }_1^{\text{amp}}>}_t$ versus breathing frequency $f_{\mathrm{p}\mathrm{B}\mathrm{ö}\mathrm{t}\mathrm{C}}$. The top panel shows the frequency computed from experimental data from head-restrained rats [Moore et al. (2013)]. The grey area highlights the range of exploratory sniffing frequencies. The bottom panel shows the frequency calculated as a function of breathing frequency for two values of $\mathrm{\Delta }{t}_{\mathrm{p}\mathrm{B}\mathrm{ö}\mathrm{t}\mathrm{C}}=20\mathrm{m}\mathrm{s}$ (bottom lines) and 70 ms (top lines), and $T_{\text{rand}}=10$ (dotted lines) and 70 ms (solid lines).

Figure 8.. Analysis of perturbations to vIRt^ret synapses and summary of brainstem control of whisking.

(A) Diagram for dynamics of the vibrissa oscillator with the vIIRt^ret under partial gephyrin degradation to weaken all ${\mathrm{G}\mathrm{A}\mathrm{B}\mathrm{A}}_{\mathrm{A}}$-ergic inhibitory inputs to vIRt^ret neurons (yellow circle). We assume that all synapses in the vIRt^ret are equally weakened by a factor of x. A. We chose $g_{\text{intra}}=12\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2,g_{\text{inter}}=20\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2,g^{\mathrm{r}\mathrm{B}}=0.5\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2$, $g_{\text{adapt}}=7\mathrm{m}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2,I_{\text{ext}}^{\mathrm{r}}=20\mu \mathrm{A}/{\mathrm{c}\mathrm{m}}^2,I_{\text{ext}}^{\mathrm{F}}=3.1\mu \mathrm{A}/{\mathrm{c}\mathrm{m}}^2,K=25$ and $N=100$ for our simulations. (B) The amplitude of the four consecutive whisks after the onset of pBötC activity, denoted by $i=1,2,3,4$, as a function of x (Equation 23). Error bars denote SD. The yellow band in B-D indicates the fraction estimated from the data in [Takatoh et al. (2022)]. (C) The times between average successive whisks ${<t_{\mathrm{w},\mathrm{i}+1}-t_{\mathrm{w},\mathrm{i}}>}_t$ as a function of as a function of x. (D) Dependence of ${⟨{\mathrm{C}\mathrm{V}}_2⟩}_i$ of vIRt^ret neurons on the value of x. (E) Summary of the circuitry that underlies the vibrissa motor plant. The vIRt drives rhythmic motion, and the shape of the waveform is set by the motor plant and feedback pathway that change the input to intrinsic and extrinsic motoneurons, and thus contact forces, upon touch. Figures updated from published summary [Bellavance et al. (2017), Kleinfeld and Deschênes (2011)].