The cerebellum shapes motions by encoding motor frequencies with precision and cross-individual uniformity

Abstract

Understanding brain behaviour encoding or designing neuroprosthetics requires identifying precise, consistent neural algorithms across individuals. However, cerebral microstructures and activities are individually variable, posing challenges for identifying precise codes. Here, despite cerebral variability, 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 confirm that deep cerebellar neurons encode frequencies using populational tuning of neuronal firing probabilities, creating cerebellar oscillations and motions with matched frequencies. The mechanism is 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 validate 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 a brain-computer interface for motor control.

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, A schematic of the vibration protocol, indicating the sequence of applied frequencies. xcorr_max, maximal cross-correlation. h, An illustration of the xcorr_max for cerebellar LFPs with compensatory motions and residual body movements (accelerometer). i, A trial-by-trial profile of the xcorr_max between cerebellar oscillations and compensatory motion (top) or accelerometer signals (bottom). j–m, The mean xcorr_max values (j and k) and Pearson correlation with Fisher’s transformation (l and m) between cerebellar oscillations and compensatory motion (j and l) or accelerometer (k and m) signals across various vibration frequencies. The statistical analysis in j was performed using the Friedman test; Friedman statistic of 26.93, P = 0.0001 (two sided). n, A three-dimensional plot of cross-correlation between cerebellar oscillations and compensatory motion, with the x axis representing the time lag, the y axis showing the vibration frequency and the z axis indicating cross-correlation values. o, A three-dimensional power spectral density (PSD) plot of the cross-correlation spectrum between cerebellar oscillations and compensatory motion, with the x axis indicating frequency, the y axis showing the vibration frequency and the z axis representing spectral power (n = 6 mice). Data are presented as mean values ± s.d. NS, not significant. **P < 0.01. Source data

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

a,b, Representative time–frequency plots (a) and PSDs (b) across various vibrating frequencies. c,d, Peak PSD amplitudes of cerebellar oscillations (Friedman statistic of 18.14, P = 0.0059 (two sided)) (c) and compensatory motions (Friedman statistic of 36, P < 0.0001 (two sided)) (d) across various vibrating frequencies. e, A 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, A linear regression analysis of the frequencies at peak PSD amplitudes for cerebellar oscillations and motor activities (36 points in 6 mice). g,h, A second-by-second linear regression analysis of the amplitude correlation (g) and frequency correlation (h) for each mouse (360 points in each mouse). i,j, The collective second-by-second analysis for all amplitudes (i) and all frequencies (j) in 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 with Fisher’s transformation and the determination coefficient (R²) of the linear regression presented in i and j, respectively (n = 6 mice, one-way ANOVA; F = 111.9, P < 0.0001). Data are presented as mean values ± s.d. *P < 0.05, **P < 0.01 and ***P < 0.001. Source data

Fig. 3. Neuronal coding for rhythmic motions.

a, A scheme of simultaneous recordings of SU neuronal activities, DCN LFPs and motion kinematics. b, A representative plot of the optetrode trajectory labelled with Dil (a representative image of one mouse; eight mice were recruited for analysis with matched cannula trajectory; Methods). c, The SU-firing rates (grey circles) and burst rates (orange circles) in DCN versus motion frequencies (n = 222 units from 8 mice). d, A scheme of the vector strength spectrum analysis. e, Vector strengths of ten SUs. 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 selection of included units. Intensity is in arbitrary units of vector strength (no unit), LFPs or motions (mV). The blue spectrum represents the mean vector strength of included units, the black spectrum represents the DCN LFP and the purple spectrum represents the motion. g, Frequency convergence of motions, LFPs and vector strength data superimposed in all trials of all mice. The top two subplots show the frequency spectrum of motion (top) and cerebellar LFP (middle). The light lines represent single trials and the heavy lines represent the averages of all trials. All peaks with sufficient prominence (Methods) detected in the vector strength spectrums throughout the expansion of the unit population (bottom). The colour gradient from green to blue reflects increasing units to calculate the vector strength spectrum. The colour depth indicates the level of prominence (n = 138 units from 8 mice). Units with a minimum spike number <10 were excluded to avoid unreliable computation of vector strength). h,i, Quantitative analysis of vector strength spectrums: the peak frequency differences to motions (h) from vector strength spectrum (left four, green to blue) or from DCN LFPs (rightmost, grey) and the SNR (i) (Friedman statistic of 213.3, P < 0.0001 (two sided)), indicating peak prominence of corresponding vector strength spectrums (n = 138 units from 8 mice). j–n, The tuning frequencies of neuronal firing probabilities via autocorrelation spectrum (j) with a representative trial (k), group analysis (l) and quantification (Δ frequency to motion (m) and SNR (n), n = 138 units from 8 mice; Friedman statistic of 292.9, P < 0.0001 (two sided)). o, A scheme of the phasic tuning of SU-firing probabilities to the instantaneous phases of motion. p,q, Representative polar plots for original (p) and shuffled (q) data. 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, Wilcoxon matched-pairs signed-rank test; W = −6,691, P < 0.0001, two sided). See Methods for detailed definitions of burst detection, vector strength and peak prominence. Data are presented as mean values ± s.d. Units represent biologically independent recordings from different neurons. ***P < 0.001. Source data

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

a, A schematic of the experimental set up 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, one-way ANOVA; F = 62.56, P < 0.0001). d–g, Vector strength analysis, including a representative example (d), group analysis (e), frequency differences between motion and vector strength spectrum peaks (f) and SNR of the spectrum peaks (g) (n = 58 units from 6 trials in 2 mice, one-way ANOVA; F = 208.2, P < 0.0001). CB, cerebellar. h, A scatter plot of the 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, Peak PSD amplitudes of cerebellar oscillations (Friedman statistic of 25.22, P = 0.0003 (two sided)) (k) and motions (Friedman statistic of 28.90, P < 0.0001 (two sided)) (l) across various stimulating frequencies. Collective data from seven trials in three 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 s intervals across all trials (2,520 points from 7 trials in 3 mice). o, Statistical analysis of the Pearson correlation with Fisher’s transformation between cerebellar LFPs and motor activity in the time domain and the determination coefficient (R²) of the linear regression presented in n (from 7 trials in 3 mice, Kruskal–Wallis test, Kruskal–Wallis statistic of 63.12, P < 0.0001, two sided). Data are presented as mean values ± s.d. Units represent biologically independent recordings from different neurons. **P < 0.01 and ***P < 0.001. Source data

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, A 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) and 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; Friedman statistic of 9.75, P = 0.0048 (two sided)). h, Statistical analysis of the correlation between cerebellar oscillation and motion in both the time domain (Pearson correlation with Fisher’s transformation) and the frequency domain (R²) (n = 8 mice; Friedman statistic of 24, P < 0.0001 (two sided)). Data are presented as mean values ± s.d. *P < 0.05 and ***P < 0.001. Source data

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

a, Optogenetic DCN stimulation with linear chirp waveform. b, A 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 (Friedman statistic of 12, P = 0.0011 (two sided)). 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; Friedman statistic of 22.2, P < 0.0001 (two sided)). 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, The activity evolution of a representative SU. i, A group analysis of correlation coefficient of DCN firings and chirp waveforms (n = 136 units from 8 mice, Wilcoxon matched-pairs signed-rank test; W = −9,018, P < 0.0001, two sided). j, Replication of the experiment with a complex chirp waveform stimulation. k, A representative time-frequency plots of the stimulation signal, cerebellar LFPs and motions. l, Frequency domain analysis (710 points in 12 mice). m, Time domain analysis (Friedman statistic of 18.17, P < 0.0001 (two-sided)). n, Statistical analysis (Friedman statistic of 32.4, P < 0.0001 (two-sided)). o, Predicted chirp points of maximal firing probability and their evolution across stimulation trials. p, The activity evolution of a representative SU. q, Group analysis of correlation coefficient of DCN firing and chirp waveforms (48 units in 12 mice; Wilcoxon matched-pairs signed-rank test; W = −1,172, P < 0.0001, two sided). Data are presented as mean values ± s.d. *P < 0.05, **P < 0.01 and ***P < 0.001. CB, cerebellum. Source data

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

a, Experimental settings of cerebellar EEG and EMG. b–e, Representative traces (b), time–frequency plots (c) and spectral diagrams of cerebellar EEG (d) and EMG (e). 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 tACS and simultaneous recording of tapping kinematics. h, The study protocol. tACS was set at the tapping frequency of 4 Hz and applied during the middle 2 min of volitional tapping. i, Frequency stability calculated from amplitude-independent kinematics (Methods). j, tACS modulation of the frequency stability of motion kinematics without a sound guide. Bidirectional modulation was observed (n = 6 subjects with 3 repeated experiments; 9 and 9 trials with increased and decreased of frequency stability, respectively). Wilcoxon signed-rank test (two sided). Increased group: P = 0.0039 (baseline versus on), 0.0547 (on versus off) and 0.6523 (baseline versus off); decreased group: P = 0.0039, 0.0742 and 0.0078, 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, Wilcoxon signed-rank test; W = 1,326, P < 0.0001, two sided). Data are presented as mean values ± s.e.m. **P < 0.01 and ***P < 0.001. Source data

Fig. 8. Summary of cerebellar frequency coding for motor kinematics.

The cerebellum encodes dynamic motor frequencies for kinematics, with supreme numerical precision and cross-individual consistency. The motor frequencies are generated by integrating neuronal firing probabilities at the populational level. The motor frequencies can be highly dynamic across time to construct non-rhythmic movements. The causality of the frequency-coding mechanism can be established by optogenetic manipulation in mice and current stimulations in humans. The cerebellar frequency codes for motion in both mice and humans are identical.