The cerebellar clock: predicting and timing somatosensory touch

Humans are adept at predicting what will happen next and when precisely it will occur. An activity as everyday as walking at a steady pace through a busy city while talking to a friend can only happen as smoothly as it does because the human brain has predicted most of the sensory feedback it will receive. It is only when the sensory feedback does not match what was expected, say, a sudden slippery spot on the pavement, that one becomes aware of the sensory feedback. The cerebellum is known to be involved in these predictions, but not much is known about the precise timing of them due to the scarcity of time-sensitive cerebellar neuroimaging studies, such as ones conducted with magnetoencephalography. We here investigated the timing of sensory expectations as they are expressed in the cerebellum using magnetoencephalography. We did this by comparing the cerebellum’s response to somatosensory omissions from regular trains of stimulation to its response to omissions from irregular trains of stimulation. This revealed that omissions following regular trains of stimulation showed higher cerebellar power in the beta band than those following irregular trains of stimulation, precisely when the omitted stimulus should have appeared. We also found evidence of cerebellar theta band activity encoding the rhythm of new sequences of stimulation Our results furthermore strongly suggest that the putamen and the thalamus mirror the cerebellum in showing higher beta band power when omissions followed regular trains of stimulation compared to when they followed irregular trains of stimulation. We interpret this as the cerebellum functioning as a clock that precisely encodes and predicts upcoming stimulation, perhaps in tandem with the putamen and thalamus. Relative to less predictable stimuli, perfectly predictable stimuli induce greater cerebellar power. This implies that the cerebellum entrains to rhythmic stimuli for the purpose of catching any deviations from that rhythm.


Abstract
Humans are adept at predicting what will happen next and when precisely it will occur. An activity as everyday as walking at a steady pace through a busy city while talking to a friend can only happen as smoothly as it does because the human brain has predicted most of the sensory feedback it will receive. It is only when the sensory feedback does not match what was expected, say, a sudden slippery spot on the pavement, that one becomes aware of the sensory feedback. The cerebellum is known to be involved in these predictions, but not much is known about the precise timing of them due to the scarcity of time-sensitive cerebellar neuroimaging studies, such as ones conducted with magnetoencephalography.
We here investigated the timing of sensory expectations as they are expressed in the cerebellum using magnetoencephalography. We did this by comparing the cerebellum's response to somatosensory omissions from regular trains of stimulation to its response to omissions from irregular trains of stimulation. This revealed that omissions following regular trains of stimulation showed higher cerebellar power in the beta band than those following irregular trains of stimulation, precisely when the omitted stimulus should have appeared. We also found evidence of cerebellar theta band activity encoding the rhythm of new sequences of stimulation Our results furthermore strongly suggest that the putamen and the thalamus mirror the cerebellum in showing higher beta band power when omissions followed regular trains of stimulation compared to when they followed irregular trains of stimulation.
We interpret this as the cerebellum functioning as a clock that precisely encodes and predicts upcoming stimulation, perhaps in tandem with the putamen and thalamus. Relative to less predictable stimuli, perfectly predictable stimuli induce greater cerebellar power. This implies that the cerebellum entrains to rhythmic stimuli for the purpose of catching any deviations from that rhythm.

Introduction
Interest in the cerebellum has surged recently. This is likely due to the realization that it does not only subserve motor behaviour and coordination as previously believed (Schmahmann, 1997). In a recent study, using functional magnetic resonance imaging (fMRI), King et al. (2019) showed that the functions of the cerebellum span many domains, such as hand movements, saccades, divided attention, verbal fluency, autobiographical recall, word comprehension, action observation, mental arithmetic, emotion processing and language processing. Thus, it is involved in many functions hitherto believed to be associated with the cerebrum. In terms of surface area, it has also recently been shown that the cerebellum is only 20% smaller than the cerebral cortex . Furthermore, it has become clear that the cerebellum actively predicts what feedback it should receive from the external world (Hull, 2020). Important questions here are what the spatiotemporal nature of these predictions is. In this study we focus on the temporal aspects, especially how the degree of irregularity in the environment affects cerebellar predictions. In a recent study (Andersen and Lundqvist, 2019), using magnetoencephalography (MEG), we found that the cerebellum was involved in maintaining and updating expectations when regularly timed somatosensory stimulation was interrupted by unexpected omissions of stimulation. We furthermore found that an evoked response in the secondary somatosensory cortex (SII) was elicited by the omitted stimuli at the exact time somatosensory stimulation was expected, despite an inter-stimulus interval (ISI) of 3 s between stimuli. This SII response has also been found with a shorter ISI, i.e. 500 ms (Naeije et al., 2018). This suggests that the temporal expectations are very precise even over long intervals. For comparison, ISIs must be shorter than at least 500 ms for auditory omissions to elicit responses (Joutsiniemi and Hari, 1989;Yabe et al., 1997).
The cerebellum has had a reputation of being impossible to study with non-invasive electrophysiological methods due to suspicions of it being too finely folded, which should result in signal cancellation, to produce signal measurable at a distance. However, we, (Andersen et al., 2020), recently reviewed the MEG and electroencephalographic (EEG) literature on cerebellar findings and found at least thirty MEG or EEG studies reported cerebellar activity. Supporting the validity of these findings, of which many did not have cerebellum as their primary target, Samuelsson et al. (2020) recently showed, using a combination of MEG signal simulation and highresolution magnetic resonance imaging (MRI) capturing the fine structure of the cerebellum, that cerebellar signals are strong enough to be detectable by state-of-the-art MEG. The time thus seems ripe for investigating the cerebellum with a targeted study. We here follow the working hypothesis that the cerebellum functions as a clock creating temporal predictions about future stimulation.
The earliest MEG evidence pointing to the cerebellum playing an important role in updating and maintaining expectations is a study by Tesche and Karhu (2000) who estimated time courses for dipolar sources that they placed in the cerebellar vermis. They found oscillatory responses to omitted stimuli in the ranges 6-12 Hz and 25-35 Hz. These were maximal around the time the stimulus should have happened, indicating a relationship to expectational processes. We, (Andersen and Lundqvist, 2019), found more direct evidence of cerebellar oscillatory responses using a dynamic imaging of coherent sources (DICS) beamformer (Gross et al., 2001). We found that the cerebellum was more strongly activated in the first repetition of a stimulus compared to the very first stimulation for the ranges 4-7 Hz and 10-30 Hz. Furthermore, we found a tonic difference between omitted stimulations, i.e. unexpected absences of stimulation, and non-stimulations, i.e. expected absences of stimulation in the range of 3-15 Hz. The tonic difference between nonstimulations and omissions was in the direction of more power in the non-stimulations in the murhythm. This may be an effect of more power in the mu-rhythm when the body is at rest (Kuhlman, 1978), as is the case for non-stimulations. Thus, to find a time-dependent cerebellar response to omissions, as is to be expected if the cerebellum functions as a clock, we reasoned that omissions of different kinds should be contrasted against one another instead of against non-stimulations.
To this end, we have created a paradigm which in many ways is similar to those of Andersen and Lundqvist (2019) Naeije et al. (2018) and Tesche and Karhu (2000). Crucially, however, we manipulated the regularity of stimulation trains by adding different levels of temporal jitter. The basic idea of the paradigm is that an underlying rhythm is introduced with three tactile stimulations following one another with the same temporal spacing between them. Then three extra stimulations follow, either jittered (irregular) or with the same temporal spacing (regular), finally followed by an omission of stimulation. The paradigm also included longer periods of non-stimulation that would follow an omission.
The main hypothesis of the current study is that the cerebellum will be more strongly activated to omissions of otherwise expected stimulation during regular than during irregular stimulation trains. We expect this difference in cerebellar activation to occur around the expected timing of the omitted stimulus following the underlying rhythm established by the first three stimuli. Based on the findings in Andersen and Lundqvist (2019) and Tesche and Karhu (2000) we expect these differences to occur in the theta and beta bands.
Following the recent findings of Andersen and Lundqvist (2019), Fardo et al. (2017), Naeije et al. (2018) and Allen et al. (2016), we also investigated other areas that are implicated in somatosensation and in updating sensory expectations, specifically the primary and secondary somatosensory cortices (SI and SII), and the inferior parietal cortex (IPC).
Finally, it must be mentioned that when speaking of timing capabilities, structures of the basal ganglia and thalamus are likely to be heavily involved (Harrington and Haaland, 1998). In fact, it has been proposed that the basal ganglia and the cerebellum together with the thalamus are nodes of an integrated network responsible for cognitive timing (Bostan and Strick, 2018;Gibbon et al., 1997), and dysfunction of this network has been implicated in Parkinson's Disease (Caligiore et al., 2016). We thus also should expect differences in the basal ganglia and the thalamus, but due to both the depth and the closed structures of neurons (Lorente de Nó, 1947) in basal ganglia structures and the thalamus, MEG is deemed anything but insensitive to these structures (Hari and Puce, 2017). However, Attal et al. (2012) (see also: (Attal and Schwartz, 2013)) argue that MEG should be at least as sensitive to the putamen of the basal ganglia as MEG is to the hippocampus, which is recently seen as a more and more viable target of MEG (Garrido et al., 2015;Ruzich et al., 2019). Attal and Schwartz (2013) also showed experimental evidence that thalamic activation may be within reach of MEG (see also : Pizzo et al., 2019). Therefore, we investigated the putamen of the basal ganglia and the thalamus.

Participants
Thirty right-handed participants volunteered to take part in the experiment (seventeen males and thirteen females, Mean Age: 26.7 y; Standard Deviation: 6.3 y; Range: 18-41 y). The experiment was approved by the local institutional review board in accordance with the Declaration of Helsinki.

Stimuli and procedure
Tactile stimulation was generated by two ring electrodes driven by an electric current generator (Stimulus Current Generator, DeMeTec GmbH, Germany). The ring electrodes were fastened to the tip of the right index finger. One was placed 10 mm below the bottom of the nail and the other 10 mm below that. Three kinds of sequences were administered. 1) a steady sequence, where all stimuli happened exactly on time; 2) a jittered sequence, where stimuli 4-6 were 5 % jittered, i.e. they happened from ~-75 ms to ~75 ms (in steps of 1 ms) relative to where the stimulus would have happened, had the sequence been steady; 3) a heavily jittered sequence, where stimuli 4-6 were 15 % jittered, i.e. they happened from ~-225 ms to ~225 ms (in steps of 1 ms) relative to where the stimulus would have happened, had the sequence been steady. For both cases, the stimulation had to be a minimum of 13 ms away from the time where the steady stimulus would have happened.
150 sequences of each kind were administered, resulting in 3,300 events of which 2,700 were stimulations and 450 were omissions. The remaining 150 were non-stimulations. Stimuli were administered in trains of six with an ISI of 1,487 ms. This made sure that stimulation would not lock to the 50 Hz power coming from electrics. After the sixth stimulus an omission of stimuli occurred and a new train of stimulation was begun. Thus, between the last stimulus of a train and the following train 2,974 ms elapsed (Fig. 1). An exception to this was that after each fifteen sequences, a train of five non-stimulations would occur, i.e. 5 x 1,487 ms of no stimulation. Thus, omissions and non-stimulations differ by whether or not an expectation was present in the time leading up to it. This structure resulted in 450 First Stimulations, 450 Repeated Stimulations (2) and 450 Repeated Stimulations (3). Note that the first three stimulations are identical in all three conditions and are thus collapsed. Furthermore, this resulted in 3 x 150, one for each of the three conditions, Repeated Stimulations (4), Repeated Stimulations (5) and Repeated Stimulations (6). Finally, it resulted in 3 x 150, one for each of the three conditions, Omitted Stimulation and a total of 150 Non-Stimulations.

Fig. 1: Stimulation sequence:
There are three conditions. A steady rhythm where all stimulations are on-time, with an inter-stimulus interval of 1,487 ms. For the two jittered conditions, the first three stimuli are on time while the last three are off-time with different amounts of jitter, randomly chosen. The omission is timed to where the stimulus would have occurred, had it been on-time, i.e. 8,922 ms after the first stimulation. Black vertical lines indicate that the stimulation is on-time, and the red lines indicate that they are jittered. The red lines furthermore indicate where the stimulus should have occurred, had it been regular.
During the stimulation procedure, participants watched a nature programme with sound from panel speakers. Participants were instructed to pay full attention to the movie and pay no attention to the stimulation of the finger.
Electro-oculography, -cardiography and myography (EOG, ECG and EMG) were recorded. For EOG, eye movements and eye blinks were recorded to monitor participants. EMG was recorded over the splenius muscles. These were used to monitor that the participant would not build up tension in the neck muscles. Lastly, for explorative purposes, respiration was measured using a respiratory belt but is not analysed here.

Preparation of participants
In preparation for the measurement, each participant had two pairs of EOG electrodes placed horizontally and vertically, respectively, around the eyes. ECG was measured by having a pair of electrodes on each collarbone. Finally, two pairs of EMG electrodes were placed on either side of the splenius muscles. Four head-position indicator (HPI) coils, two behind the ears and two on the forehead, were placed on the participants. Subsequently, each participant had their head shape digitized using a Polhemus FASTRAK. Three fiducial points, the nasion and the left and right preauricular points, were digitized along with the positions of the HPI-coils. Furthermore, around 200 extra points digitizing the head shape of each participant were acquired. Participants were subsequently placed in the supine position of the MEG system and great care was taken, such that they would lie comfortably in the scanner, thus preventing neck tension.

Acquisition of data
Data was recorded on an Elekta Neuromag TRIUX system inside a magnetically shielded room (Vacuumschmelze Ak3b) at a sampling frequency of 1,000 Hz. As data was acquired online, lowpass and high-pass filters were applied, at 330 Hz and 0.1 Hz respectively.

Processing of MEG data
We analysed evoked responses and oscillatory responses. MaxFilter was not applied to preserve the rank of the data. The four first signal space projections (SSPs) were projected out for the sensor space analyses due to the strong influence of gradients on the magnetometers.
For the evoked responses, we low-pass filtered the data at 40 Hz (finite impulse response; zerophase; -6 dB cutoff frequency: 45 Hz; order: 198) and then cut the raw data into epochs of 800 ms, 200 ms pre-stimulus and 600 ms post-stimulus. The epochs were demeaned using the pre-stimulus period. Segments of data including magnetometer responses greater than 4 pT or gradiometer responses greater than 400 pT/m. Subsequently, epochs were averaged to create evoked responses for each of the trial types.
For the oscillatory responses, the data were band-pass filtered into the theta and beta bands, both zero-phase with a -6 dB cutoff frequency: from 4-7 Hz (Theta) and 14-30 Hz (Beta). A Hilbert transform was applied to both bands and the data were cut into epochs of 1,500 ms, 400 ms prestimulus and 1,100 ms post-stimulus for the stimulations and 750 ms pre-stimulus and 750 ms poststimulus for the omissions. The stimulation epochs were demeaned using the pre-stimulus interval up to -50 ms, whereas the omission epochs were demeaned using the whole time-interval. The difference in demeaning was due to the lack of evoked responses for omission epochs. Segments of data including magnetometer responses greater than 4 pT or gradiometer responses greater than 400 pT/m were rejected. Epochs were not rejected based on EOG and EMG due to beamformers (see below) being good at suppressing the artefacts arising from eye blinks and movements. We calculated the envelopes of the Hilbert transformed data, and, finally, to minimise the effect of outliers on the average of epochs, for each contrast of interest, we converted that contrast to a zscore, using a Wilcoxon signed-rank test entering each condition of the contrast into the test.

Source reconstruction
For both the evoked responses and the oscillatory responses, a linearly constrained minimum variance (LCMV) beamformer (van Veen et al., 1997) was applied. For the former, it was applied to evoked responses and for the latter to the Hilbert transformed epochs.
We acquired sagittal T1 weighted 3D images for each participant using a Siemens Magnetom Prisma 3T MRI. The pulse sequence parameters were: 1 mm isotropic resolution; field of view: 256 mm x 256 mm; 192 slices; slice thickness: 1 mm; bandwidth per pixel: 290 Hz/pixel; flip angle: 9°, inversion time (TI): 1,100 ms; echo time (TE): 2.61 ms; repetition time (TR): 2,300 ms. Based on these images, we did a full segmentation of the head and the brain using FreeSurfer Fischl et al., 1999). Subsequently we delineated the head and brain surfaces using the watershed algorithm from MNE-C (Gramfort et al., 2013). We created a volumetric source space from the brain surface with sources 7.5 mm apart from one another (~4000 sources). Single compartment boundary element method (BEM) solutions (volume conduction models) were calculated from the brain surfaces. For each participant, the T1 was registered to the participant's head shape with the fiducials and head shape points acquired with Polhemus FASTRAK. Finally, with the co-registered T1, we created two forward models based on the volume conduction model, the positions of MEG sensors, and the volumetric source space.
For estimating source time courses, the data covariance matrix was estimated based on the poststimulus period (from 0 ms to 600 ms) for the evoked responses. For the oscillatory responses, the data covariance matrix was estimated based on the post-stimulus time period (from 0 ms to 1,100 ms) for the stimulations, and for the omissions the data covariance matrix was estimated based on the whole time period (-750 ms to 750 ms). No regularisation was applied when estimating the filter weights. For the filters, we chose the source orientation that would maximise source output. Finally, for the oscillatory responses, we took the absolute value of the complex-valued LCMV-source-timecourses, and averaged over all the epochs to acquire an averaged source time course for both bands (4-7 Hz and 14-30 Hz). For the evoked responses, we also used the absolute value of the estimated source time courses.
For all LCMVs only magnetometers were used since these are the sensors most sensitive to deep sources. These source reconstructions were morphed onto fsaverage, a common template from FreeSurfer. Note, that signal space projection vectors were not applied on the data to be beamformed since beamformers suppress low-rank external noise well (Sekihara and Nagarajan, 2008).
To investigate contrasts in an unbiased manner, the spatial filter of the LCMV beamformer was estimated based on the covariance of the combined data. Subsequently, this filter was used to estimate the power for each condition in the contrast and finally the contrast for the oscillatory responses was calculated as: (condition_1 -condition_2) / (condition_1 + condition_2), meaning that the contrast expresses the difference in power as a percentage of total power in the two conditions. For the evoked responses, we looked at the difference (condition_1 -condition_2).

Statistical analysis
For the oscillatory responses, we focused on the comparisons between the Omissions, with the difference between Omission 0 and Omission 15 expected to be the strongest. Also the comparison between First Stimulation and Repeated Stimulation (2) was interesting due to the difference in expectation and our earlier findings of a difference in the theta band (Andersen and Lundqvist, 2019). We used a cluster permutation test (Maris and Oostenveld, 2007) to address the multiple comparisons problem arising from doing all-sensor and whole-brain analyses. The analyses done in sensor space were done on the whole time range (-400 ms to 1,100 ms for the stimulations and -750 ms to 750 ms for the omission). We used this to inform the time range of our subsequent tests for the whole-brain analysis. The null hypothesis of such a test is that data are exchangeable, meaning colloquially speaking that our labelling does not matter. The alternative hypothesis then is that the labelling does matter. The procedure was as follows. First we conducted a t-test for each of the time-sensor or time-source combinations (e.g. for the beamformer: 801 time samples and 4,342 sources; > 3 million tests). Then spatio-temporal clusters were formed based on connecting neighbouring sensors/sources and neighbouring time points. Only time-sensor/source combinations that were significant at α = 0.05 were included in clusters. 1024 permutations were then run where the condition labels, e.g. Omission_0 and Omission_15, were shuffled between trials. Using the same procedure as for the correctly labelled trials, spatio-temporal clusters were formed from the significant tests for each of the permutations. A null hypothesis distribution was sampled by adding the value of the largest cluster to the distribution for each permutation. The largest cluster in each permutation was defined as the one that had the largest sum of absolute t-values. For each cluster in the correctly labelled trials, the likelihood of that cluster, or one more extreme, was found, its pvalue, given the sampled null hypothesis distribution. The null hypothesis that the data is exchangeable was then rejected if the largest cluster in the correctly labelled data was associated with a p-value that is lower than the alpha value, which we set at 0.05.
For the evoked responses, we did the cluster permutation tests on the source time courses directly, since the SI and SII responses were expected. We focused on the comparisons between First Stimulation and Repeated Stimulation (2) and between Omitted Stimulation and Non-Stimulation following Andersen and Lundqvist (2019). We tested on the full time range (-200 ms to 600 ms).

Evoked responses
The electric stimulations gave rise to the expected SI and SII activations ( Surprisingly, no omission responses were found ( Fig. 2A). This is in contrast to two recent studies (Andersen and Lundqvist, 2019;Naeije et al., 2018), in which evoked responses source localized to SII were found peaking around 130 ms. A difference between these two studies and the B A C current, is that we in the current study used electric stimulation, whereas the two other studies used an inflatable membrane that pressed on the finger when inflated. Differences were found between First Stimulation and Repeated Stimulation (2) for the source time courses when using the cluster permutation procedure described above in the time range -200 ms to 600 ms (p BIGGEST_CLUSTER = 0.00098). No differences were found between Omission 0 and Non-Stimulation (p BIGGEST_CLUSTER = 0.73).

Oscillatory responses
For the oscillatory responses, we focused on the theta and beta band responses, following the results of Andersen and Lundqvist (2019) and Tesche and Karhu (2000). For the stimulations, the main contrast of interest was between First Stimulation and Repeated Stimulation (2). In the beta band, the null hypothesis could not be rejected (p BIGGEST_CLUSTER = 0.72). In the theta band, we found that 47 magnetometers formed part of two spatiotemporal clusters with p-values < 0.05 (p BIGGEST_CLUSTER = 0.00098) (Fig. 3A). Given that the responses peaked around 200 ms, we centred the subsequent source cluster analysis on the window of -200 ms to 600 ms. We found one cluster with a p-value < 0.05 (p BIGGEST_CLUSTER = 0.00098). The maximal value was found for the left inferior parietal cortex at 145 ms (Fig. 3D), beginning at 0 ms (Fig. 3F). A left cerebellar activation was also found after 290 ms (Fig. 3C),. For the comparison between Repeated Stimulation (2) and Repeated Stimulation (3) the null hypothesis could not be rejected (p BIGGEST_CLUSTER = 0.12). For descriptive purposes, we reconstructed all the first three stimulations based on a common spatial filter and estimated time courses for the sources shown in Fig. 3CD (Fig. 3EF). These show the cerebellar and inferior parietal cortical activity for each of the conditions in separation. For the omissions, the main contrast of interest was between Omission 0 and Omission 15, as we expected the cerebellar contrast to be the strongest there. In the theta band, the null hypothesis could not be rejected (p BIGGEST_CLUSTER = 0.25). In the beta band, we found that 14 magnetometers formed part of three spatiotemporal clusters with p-values < 0.05 (p BIGGEST_CLUSTER = 0.014 ) (Fig.  4A). The sensor topography (Fig. 4B) was compatible with an underlying cerebellar source on the right side. For the source reconstruction, we found one cluster with a p-value < 0.05 (p BIGGEST_CLUSTER = 0.0078) testing on the time interval from -400 ms to 400 ms. A right cerebellar source was part of the cluster (Fig. 4C), and the maximally responding source 0 ms was in the putamen (Fig. 4D). Running similar tests for Omission 0 vs Omission 5 and Omission 5 vs Omission 15 resulted, however, in not being able to reject the null hypothesis (p BIGGEST_CLUSTER = 0.18 and p BIGGEST_CLUSTER = 0.11 respectively). Again, for descriptive purposes, we reconstructed all three types of omissions based on a common spatial filter and estimated time courses for the sources shown in Fig. 4CD  (Fig. 4EF). These show that cerebellar and putamen activity starts increasing around -250 ms of stimulation when the omission is preceded by a regular train of stimulation (Omission 0) and starts decreasing again around 0 ms. The opposite pattern emerges for the omissions preceded by irregular trains of stimulation (Omission 5 and Omission 15), which decrease from -250 ms and start increasing from 0 ms. We investigated the other areas of interest besides cerebellum and putamen. For the contrast between Omission 0 and Omission 15 in the beta band, we also found differences in the thalamus, the inferior parietal cortex, SI and SII (Fig. 5) In summary, we found that non-cortical structures, i.e. cerebellum, putamen and thalamus, followed the expected timing, i.e. the contrast peaking at ~0 ms, whereas cortical structures such as SI and inferior parietal cortex peaked before the expected stimulus (~-250 ms) and SII peaking after the expected stimulus (~250 ms).

for the individual time courses.
We then investigated the same areas for the contrast between First Stimulation and Repeated Stimulation (2) for the theta band (Fig. 6). Here, we found differences in the very same areas. Except for the cerebellum, all differences were expressed around the peak of the left inferior parietal cortex (~145 ms) (Fig. 3D&F&H). The cerebellar differences peaked around 285 ms for the left side and around 500 ms for the right side.  (Tzourio-Mazoyer et al., 2002) with the exception of SII, which is from the Harvard Oxford cortical atlas (Desikan et al., 2006). See Supplementary Fig. 2

Discussion
In this study we investigated if and how the cerebellum is involved in timing expectations. We presented trains of electrical stimulation that were either regularly timed or irregularly timed and followed by an omission of stimulation. We found that our hypothesis, that cerebellum more strongly activates to omissions appearing in a regularly timed train than in an irregularly timed train, to be corroborated by our beta band findings, showing that the contrast and the omission activity related to the regular train was at their strongest close to 0 ms (Fig. 4E & 4G).

A cerebellar clock (the beta band)
That the differences in cerebellum is maximal around the expected time of stimulation indicates that cerebellar power reflects the strength of the temporal expectation, which theoretically should be strongest in the predictable condition. In this regard, the cerebellum functions as a clock keeping track of when sensory information is supposed to arrive. On the contrary, when there is temporal uncertainty regarding when the stimulus should arrive cerebellar activity decreases around the time of the expected stimulation, and for the most uncertain one (Omission 15) most so (Fig. 4E). Cerebellar power thus reflects the prediction of an upcoming somatosensory stimulation, as beta band power is expected to increase when predicted information is expected to arrive (Engel and Fries, 2010) (Fig. 4E). The cerebellar findings here are also similar to findings from the auditory domain, where it has been shown that beta band power in the auditory cortex synchronizes with the rhythm of auditory stimuli, peaking at their occurrence and then decreasing again (Arnal and Giraud, 2012;Fujioka et al., 2012).
Interestingly, the most irregular condition (Omission 15) showed a dip for oscillatory power around the time, the stimulation should have occurred, had the pattern been regular (Fig. 4G). In this condition, a stimulation is expected, but its exact timing is uncertain. Given that a function of prediction is to inhibit the processing of the expected stimulus, we can interpret this as the cerebellum uninhibiting the processing of the spatially expected, but temporally unexpected, stimulus. Interestingly, the dip begins just before the earliest time that the stimulation could occur (~-225 ms) and is back to pre-expected-stimulus levels at around 225 ms. This suggests that the clock is even tracking the range within which the stimulus should occur. This possibility is seemingly contradicted by the observation that the slightly jittered condition (Fig. 1) also starts decreasing around 225 ms (5 % of 1,487 ms = 74 ms). This can however be explained by us having used an interleaved design; i.e. when the cerebellum detects jitter in the stimulation, there is no principled way of deciding, in the current design, whether a given stimulation train is 5 % jittered because any 5 % jittered train is compatible with being a 15 % jittered train. To test whether the decrease adapts this precisely to the temporal uncertainty, a blocked design could prove useful in subsequent studies.
Our results furthermore indicated that the putamen and thalamus are more strongly activated for the omission preceded by regular trains than omissions preceded by irregular trains. This fits well with the knowledge that the thalamus and basal ganglia are involved in timing as a recent metaanalysis of fMRI timing studies showed (Teghil et al., 2019). Of notice is that they follow the timing of the cerebellum, peaking around 0 ms, indicative of a network of cerebellum, putamen and thalamus (Bostan and Strick, 2018;Caligiore et al., 2016;Gibbon et al., 1997). For these, we also saw the dip in power for the most irregular condition (Omission 15) (Fig. 4H, Supplementary Fig.  1).
Importantly though, activations of putamen and thalamus found using MEG must be treated with caution. Both are not optimal targets for MEG, due to their deep locations and the closed field of their neurons (Lorente de Nó, 1947). This results in a weak signal when measured outside the head that makes it challenging to say with confidence which deep source produced the MEG signal. However, theoretical work and simulation studies show that MEG should be sensitive to thalamus and putamen (Attal et al., 2012;Attal and Schwartz, 2013), albeit at a degree much lower than its sensitivity to the neocortex. In terms of experimental evidence, it was recently shown that MEG can retrieve patterns of activation from the putamen and thalamus (Pizzo et al., 2019). This theoretical and experimental evidence together with the knowledge that the present paradigm should elicit activity in the basal ganglia and the thalamus combined with the high number of omission trials (N=150) mean that we believe that the present evidence is supportive, but not conclusive, of actual putamen and thalamus activations measured by MEG.
In contrast to the cerebellum, the SI and inferior parietal cortical contrasts peaked before (Fig. 5) the expected arrival of the stimulus and the SII contrast peaked after the expected arrival of the stimulus. The SI and inferior parietal cortex differences may indicate a dampening of cortical activity before the onset of an expected stimulus, compatible with the decrease in theta band activity for repeated stimulations (Figs. 3 & 6). The SII activity may reflect an update of the current state of affairs, i.e. that the train of stimulation has been broken. Similar results have been reported for omissions in the auditory domain, with beta band power increasing in the time interval after the expected, but omitted, stimulus in the auditory cortex (Fujioka et al., 2009).
We unexpectedly did not find an SII evoked response to the omissions (Fig. 2), despite two recent studies having reported this (Andersen and Lundqvist, 2019;Naeije et al., 2018). Both these studies used tactile stimulation by inflating a membrane fastened to the fingertips of participants and found an SII response around 150 ms time-locked to the expected, but omitted, stimulation. A possible explanation for not finding such a response in the current study is that the instantaneousness of the electrical stimulation used in the current study relative to the longer extended touch of the membrane makes it harder for the brain to time-lock to the exact time point. This would have the consequence that evoked analyses would not be sensitive to the SII response, but that the Hilbert transformation strategy that we applied would, given its sensitivity to responses that are not precisely time-locked. Thus, it is possible that the reported SII response (Fig. 5) is similar to the response reported by Andersen and Lundqvist (2019) and by Naeije et al. (2018).
To summarize, the current findings provide evidence that the cerebellum, together with the putamen and thalamus, tracks the timing of upcoming stimulation, as reflected by beta band power. This fits the interpretation of the beta band as predicting the when of upcoming stimulation (Arnal and Giraud, 2012), and supports the notion of the cerebellum as a clock that keeps track of upcoming stimulation. Cortical structures such as the SI, SII and the inferior parietal cortex are also timed relative to the upcoming stimulation. However, the SI and the inferior parietal cortex seem to be downregulating before the expected stimulus, and SII activates after the omitted stimulus, possibly reflecting an updating of the current state of affairs. Notably, the omission-related SII difference is strongest in the right hemisphere, mirroring the evoked difference found by Andersen and Lundqvist (2019).

Encoding of expectations (theta band)
We furthermore found in the theta band that the cerebellum reacted more strongly to the first stimulation of a train compared to the subsequent stimulation around 280 ms after stimulation onset. This may reflect the unexpected nature of the first stimulation compared to the subsequent stimulation. This is similar to what has been found using electrical stimulation (Tesche and Karhu, 2000; see also : Naeije et al., 2018), but opposite to what has been found using tactile stimulation (Andersen and Lundqvist, 2019). Future studies will have to tell whether this difference is dependent on the type of stimulation. Similarly to the omission contrast, cortical areas, i.e. SI, SII, and inferior parietal cortex, and non-cortical areas, i.e., cerebellum, putamen and thalamus, were found to differ in activation between the first and the subsequent stimulation (Fig. 6). Except for the cerebellum, all contrasts peaked around 145 ms (Fig. 3F), whereas cerebellum peaked later, around 285 ms (Fig. 3E). All the cortical areas were expressed most strongly on the contralateral side of the stimulation as expected. The difference found for the SII is likely to be the same activity caught by the evoked analysis (Fig. 1C), as the evoked response is in the theta range. The later cerebellar activations are potentially an updating of the new state of affairs (Engel and Fries, 2010), i.e., that the stimulation train has been interrupted.

Cerebellum's role in forming expectations
The current study reveals that cerebellum (along with its potential network members, thalamus and putamen (Bostan and Strick, 2018)) was the only area that peaked around 0 ms when comparing omissions (Fig. 5). We interpret this as the cerebellum clocking that stimuli arrive as expected, possibly in unison with the putamen and the thalamus. This clocking activity is expressed in the beta band. Similar clocking activity has been revealed in the auditory domain for the beta band (Arnal and Giraud, 2012;Fujioka et al., 2012;Ruiz et al., 2017).
We also found cerebellar activity in the theta band when comparing the first stimulation to the subsequent stimulation. In contrast to the beta band activity, the theta cerebellar activity was unique in peaking later (~285 ms) than the first stimulation related activations (~145 ms) (Fig. 6). We interpret this as the cerebellum updating its clock for subsequent stimulation. Of note is also that, in contrast to the omission-related beta band activity, the putamen and the thalamus peak out of sync with the cerebellum. Our findings are similar to recent findings (Dave et al., 2020), where transcranial magnetic stimulation was used to knock out cerebellar theta and cerebellar beta respectively. Cerebellar theta was found to be related to encoding episodic memories and cerebellar beta was found to be related to semantic predictions related to those memories. We can thus interpret our current findings as cerebellar theta setting the clock and cerebellar beta checking whether the clock is set correctly.

Control analyses
An alternative explanation for why we found differences in cerebellar activity when comparing omissions following regular and irregular trains of stimulation may be that the proposed cerebellar clock is based on local timelocking rather than global timelocking. By global timelocking, we mean that the timing of the expected, but omitted, stimulus is set to 8,922 ms (six times the ISI (1,487 ms)) after the first stimulation of the train (Fig. 1). By local timelocking, we mean that the timing of expected, but omitted, stimulus is set to the temporal distance between the fifth and sixth stimulation, whatever that may be in each particular stimulation train. Using local timelocking, we were also able to reject the null hypothesis (p BIGGEST_CLUSTER = 0.0029), and we found similar activations of the cerebellum and the putamen (Supplementary Fig. 3). We also used mean timelocking, where we calculated mean distance between the last three stimulations (Fig. 1). Using mean timelocking, we found similar activations of the cerebellum and the putamen ( Supplementary  Fig. 4) (p BIGGEST_CLUSTER = 0.0049). The way we timelocked the data in the main analysis (Fig. 4) thus does not change the results significantly.
It is also possible that the contributions of the gradient could not be filtered by the beamformer. We therefore re-ran the beamformer analysis on data where the signal space projections had been applied. This again resulted in very similar results (Omission 0 vs. Omission 15; p BIGGEST_CLUSTER = 0.0049; Supplementary Fig. 5).

Conclusion
In conclusion we find that the cerebellum predicts the timing of prospective somatosensory stimuli, functioning like a clock. This is evidenced by cerebellar beta-band power increasing and subsequently decreasing around the expected timing of stimulation, when stimuli are perfectly predictable, whereas when stimuli are less than perfectly predictable, cerebellar beta-band power decreases around around the expected time of stimulation, indicating that cerebellar beta band activity is related to prediction. This decrease is suggestive evidence of the clock tracking the uncertainty range of the expected upcoming stimulation, but further studies are needed for ascertaining that the range is also encoded. Also of note, we found evidence of cerebellar theta-band activity potentially encoding memories on which the subsequent predictions can be based. We interpret this as cerebellar theta activity reflecting setting the proposed cerebellar clock and the cerebellar beta activity reflecting checking the prediction of the clock. Furthermore, we find intriguing evidence of putamen and thalamus activation, tracking the timing of somatosensory stimuli as well. This fits well with knowledge from other domains, such as functional magnetic resonance imaging. However, given the low sensitivity of magnetoencephalography to the putamen and the thalamus, further research is warranted.