The cerebellum computes frequency dynamics for motions with numerical precision and cross-individual uniformity

Abstract

Cross-individual variability is considered the essence of biology, preventing precise mathematical descriptions of biological motion^1-7 like the physics law of motion. Here we report that the cerebellum shapes motor kinematics by encoding dynamic motor frequencies with remarkable numerical precision and cross-individual uniformity. Using in-vivo electrophysiology and optogenetics in mice, we confirmed that deep cerebellar neurons encoded frequencies via populational tuning of neuronal firing probabilities, creating cerebellar oscillations and motions with matched frequencies. The mechanism was consistently presented in self-generated rhythmic and non-rhythmic motions triggered by a vibrational platform, or skilled tongue movements of licking in all tested mice with cross-individual uniformity. The precision and uniformity allowed us to engineer complex motor kinematics with designed frequencies. We further validated the frequency-coding function of the human cerebellum using cerebellar electroencephalography recordings and alternating-current stimulation during voluntary tapping tasks. Our findings reveal a cerebellar algorithm for motor kinematics with precision and uniformity, the mathematical foundation for brain-computer interface for motor control.

Keywords: Cerebellum; Electroencephalogram; Frequency; Motor control; Motor kinematics; Neuronal coding; Oscillations.

Fig. 1 |. Self-generated cerebellar oscillations in compensatory motions.

(a) An experimental setting of the vibration platform generating horizontal sinusoidal motions. (b) Representative traces for active compensatory motion, calculated as signals of a head-mounted accelerometer minus the platform vibrations. (c–f) Representative time-frequency plots of vibrations (c), head-mounted accelerometer signals (d), compensatory motions (e), and cerebellar oscillations (f) during 16-Hz vibrations. (g) Schematic of the vibration protocol, indicating the sequence of applied frequencies. (h) Illustration of maximal cross-correlation for cerebellar LFPs with compensatory motions and residual body movements (accelerometer). (i-k) Trial-by-trial (i) and group analysis (j–k) of cross-correlation from cerebellar oscillations versus compensatory motions or accelerometer signals at various vibrating frequencies. Cerebellar oscillations exhibited significant correlations with compensatory motions (j) but not with residual body movements (k) (n = 6 mice). Error bars denote S.D. **p < 0.01, One-way ANOVA.

Fig. 2 |. Correlation of cerebellar oscillations and rhythmic motions in the frequency domain.

(a-b) Representative time-frequency plots (a) and power spectral density (PSD) (b) across various vibrating frequencies. (c-d) Peak PSD amplitudes of cerebellar oscillations (c) and compensatory motions (d) across various vibrating frequencies. (e) Linear regression analysis of peak PSD amplitudes between cerebellar oscillations and compensatory motor movements. The solid red line represents the best-fit linear model, while the dashed red lines indicate the 95% confidence bounds (36 points in 6 mice). (f) Linear regression analysis of the frequencies at peak PSD amplitudes for cerebellar oscillations and motor activities (36 points in 6 mice). (g-h) Second-by-second linear regression analysis for each mouse (360 points in each mouse). (i-j) Collective second-by-second analysis for all mice combined (2,160 points in 6 mice). (k) Statistical analysis of the correlation between cerebellar LFPs and motor activity in both the time domain and frequency domain, using Pearson correlation and the determination coefficient (R²) of the linear regression presented in Fig. 2i–j respectively (n = 6 mice). Error bars denote S.D. p* <0.05, **p < 0.01, One-way ANOVA.

Fig. 3 |. Neuronal coding for rhythmic motions.

(a) Scheme of simultaneous recordings of single-unit (SU) neuronal activities, DCN LFPs, and motion kinematics. (b) A representative plot of the optetrode trajectory labeled with Dil (see Methods). (c) SU-firing rates (gray circles) and burst rates (orange circles) in DCN versus motion frequencies (n = 222 units from 8 mice). (d) Scheme of the vector strength spectrum analysis. (e) Vector strengths of 10 single units. (f) Frequency convergence of the vector strength of a representative trial during 16-Hz vibration. The vector strength spectrum peaks converged to the motion frequency throughout the random recruitment of units. Intensity is in arbitrary units of vector strength (no unit), LFPs, or motions (mV). The blue spectrum represents the mean vector strength of recruited units, the black spectrum represents the DCN LFP, and the purple spectrum represents the motion. (g) Frequency convergence of motions, LFPs, and vector strengths in all trials. The top two subplots showed the frequency spectrum of motion (top) and cerebellar LFP (middle). Light lines represent single trials, and heavy lines represent the averages of all trials. The bottom figure showed all peaks with sufficient prominence (see methods) detected in the vector strength spectrums throughout the random recruitment of units. The color gradient from green to blue reflected increasing units recruited to calculate the vector strength spectrum. The color depth indicated the level of prominence (n= 138 units from 8 mice. Units with minimum spike number < 10 were excluded to avoid unreliable computation of vector strength). (h-i) Quantitative analysis of vector strength spectrums. Peak frequency differences to motions (h) from vector strength spectrum (left four, green to blue) or from DCN LFPs (rightmost, gray), and the signal-to-noise ratio (SNR, Fig. 3i), indicating peak significance of corresponding vector strength spectrums. (j-n) The tuning frequencies of neuronal firing probabilities via autocorrelation spectrum (j) with a representative trial (k), group analysis (l), and quantification (m-n). (o) Scheme of the phasic tuning of SU firing probabilities to the instantaneous phases of motion. (p-q) Representative polar plots. DCN neurons had a greater phasic bias to the phase of motion, quantified by the polarity index. (r-s) Group analysis of cumulative probabilities (r) and values (s) of polarity indexes. DCN neurons revealed stronger phasic tuning to 16-Hz compensatory motion at the populational level (n = 138 units from 8 mice). See methods for detailed definitions of burst detection, vector strength, and peak prominence. Error bars denote S.D. ***p < 0.001, One-way ANOVA (i, n), Wilcoxon matched-pairs signed rank test (s).

Fig. 4 |. Cerebellar and motor responses to optogenetic DCN stimulation at multiple frequencies.

(A) Schematic of the experimental setup and representative histology of channelrhodopsin-2 (ChR2)-expressing DCN. (B) Representative traces showing SU firing rates (top) and their modulation during 16-Hz optogenetic stimulation of the DCN (bottom). (C) Statistical analysis of SU firing rates across different phases of the 16-Hz stimulation cycle (n = 58 units from 6 trials in 2 mice). (D–G) Vector strength analysis, including a representative example (D), group analysis (E), frequency differences between motion and vector strength spectrum peaks (F), and signal-to-noise ratio of the spectrum peaks (G). (H) The scatter plot of peak cerebellar LFP frequencies against combined vector strength spectrum peaks under various stimulating frequencies. (I-J) Representative time-frequency plots (I) and spectral diagrams (J) of optogenetically driven cerebellar oscillations and corresponding motor activities. (K-M) Collective data from 7 trials in 3 mice showing the close correspondence between cerebellar oscillatory and motor frequencies (M). (N) Scatter plots of the amplitudes (left) and frequencies (right) of cerebellar LFPs and motor activity, compiled from 1-second intervals across all trials (2,520 points from 7 trials in 3 mice). (O) Statistical analysis of the correlation between cerebellar LFPs and motor activity in the time domain and the determination coefficient (R²) of the linear regression presented in Fig. 4N. Error bar denotes S.D. ***p<0.001. Kruskal-Wallis test and one-way ANOVA.

Fig. 5 |. Non-rhythmic cerebellar oscillations and motor kinematics induced by linear chirp vibrations.

(a) The experimental settings and platform vibrations with constantly changing chirp waveform. (b) Schematic representation of vibration protocol and the time-frequency plot of the vibration signals. (c) Representative traces for compensatory motions. (d-e) Frequency domain analysis. A representative time-frequency plot of cerebellar LFPs, motions, and accelerometer signals (ACC) (d). Linear regression analysis of second-by-second amplitudes and frequencies between the cerebellar LFPs and motions (e, 2,400 points from 80 trials in 8 mice). (f-g) Time domain analysis. Trial-by-trial (f) and group analysis (g) of cross-correlation for cerebellar oscillations between compensatory motions and residual body movements (accelerometer). (h) Statistical analysis of the correlation between cerebellar oscillation and motion in both the time domain (Pearson correlation) and the frequency domain (R²) (n = 8 mice). Error bars denote S.D. *p < 0.05, ***p < 0.001, One-way ANOVA.

Fig. 6 |. Non-rhythmic cerebellar oscillations and motor kinematics induced by optogenetic stimulation.

(a) Optogenetic DCN stimulation with linear chirp waveform. (b) Representative time-frequency plot of stimulating signals, cerebellar LFPs, and motions. (c) Frequency domain analysis, linear regression analysis of second-by-second amplitudes and frequencies between the cerebellar LFPs and motions (239 points in 8 mice). (d) Time domain analysis. Trial-by-trial (left) and group analysis (right) of cross-correlation for cerebellar LFPs between motions. (e) Statistical analysis of the correlation between cerebellar oscillations and motions in both the time domain (Pearson correlation) and frequency domain (R²) (n=8 mice). (f) SU activities of DCN with linear chirp-wave stimulation. (g) Predicted chirp points of maximal firing probability and their evolution across stimulation trials (defined by the number of peaks of chirp waves). (h) Activity evolution of a representative SU. (i) Group analysis of correlation coefficient of DCN firings and chirp waveforms (n = 136 units from 8 mice). (j-q) Complex chirp waveform stimulation (l, 710 points in 12 mice; q, 48 units in 12 mice). Error bars denote S.D. *p < 0.05, **p<0.01, ***p < 0.001, One-way ANOVA (D, E, M, N), Wilcoxon matched-pairs signed rank test (i, q).

Fig. 7 |. Cerebellar oscillations and their frequency modulation during volitional tapping of health subjects.

(a) Experimental settings of cerebellar EEG and electromyography (EMG). (b-e) Representative traces (b), time-frequency plots (c) and, spectral diagram (d-e) of cerebellar and EMG. (f) Linear regression analysis of second-by-second amplitudes and frequencies of cerebellar oscillations and EMG activities at the tapping frequencies (1,286 points, n = 10 subjects). (g) Cerebellar transcranial alternating current stimulation (tACS) and simultaneous recording of tapping kinematics. (h) Study protocol. tACS was set at the tapping frequency of 4 Hz and applied during the middle 2 minutes of volitional tapping. (i) Frequency stability calculated from amplitude-independent kinematics (see methods). (j) tACS modulation of the frequency stability of motion kinematics without a sound guide. Bi-directional modulation was observed (n = 6 subjects with 3 repeated experiments; 9 and 9 trials with increased and decreased of frequency stability, respectively). (k) tACS modulation of the frequency stability of motion kinematics with a sound guide. No significant modulation was observed. (l) Cross-correlation (xCorr) peaks between tapping kinematics and tACS waveform. Values in the sound-on period were significantly higher than the sound-off period (the same 18 trials in 6 subjects). Error bars denote S.E.M.. **p < 0.01, ***p < 0.001, Wilcoxon signed-rank test.