Abstract
fMRI research has revealed that cerebellar lobule VIIb/VIIIa exhibits load-dependent activity that increases with the number of items held in visual working memory (VWM). However, it remains unclear whether these cerebellar responses reflect processes specific to VWM or more general visual attentional mechanisms. To investigate this question, we examined whether cerebellar activity during the delay period of a VWM task is selective for stimuli held in working memory. A sample of male and female human subjects performed a VWM continuous report task in which they were retroactively cued to remember the direction of motion of moving dot stimuli. Cerebellar lobule VIIb/VIIIa delay-period activation accurately decoded the direction of the remembered stimulus, as did frontal and parietal regions of the dorsal attention network. Arguing against a motor explanation, no other cerebellar area exhibited stimulus specificity, including the oculomotor vermis, a key area associated with eye movement control. Finer-scale analysis revealed that the medial portion of lobule VIIb and to a lesser degree the lateral most portion of lobules VIIb and VIIIa, which exhibit robust resting state connectivity with frontal and parietal regions of the dorsal attention network, encoded the identity of the remembered stimulus, while intermediate portions of lobule VIIb/VIIIa did not. These findings of stimulus-specific coding of VWM within lobule VIIb/VIIIa indicate for the first time that the distributed network responsible for the encoding and maintenance of mnemonic representations extends to the cerebellum.
SIGNIFICANCE STATEMENT There is considerable debate concerning where in the brain the contents of visual working memory (VWM) are stored. To date, this literature has primarily focused on the role of regions located within cerebral cortex. There is growing evidence for cerebellar involvement in higher-order cognitive functions including working memory. While the cerebellum has been previously shown to be recruited by VWM paradigms, it is unclear whether any portion of cerebellum actively encodes and maintains mnemonic representations. The present study demonstrates that cerebellar lobule VIIb/VIIIa activity patterns are selective for remembered stimuli and that this selectivity persists in the absence of perceptual input. These findings provide novel evidence for the participation of cerebellar structures in the persistent storage of visual information.
Introduction
Visual working memory (VWM) can be broadly defined as the temporary maintenance of visual information independent of the constant influx of perceptual input. Investigating the neural mechanisms of VWM is of considerable interest to the neuroscience community as performance on VWM tasks has been shown to be highly correlated with broader measures of cognitive ability such as fluid intelligence, reading comprehension, and scholastic aptitude (Engle et al., 1999; Engle, 2002; Fukuda et al., 2010). The question of where working memory contents are stored in the brain is the subject of ongoing debate. Based on electrophysiological recordings in non-human primates and neuroimaging in humans, it has long been asserted that parietal and prefrontal cortices support working memory maintenance (Fuster and Alexander, 1971; Funahashi et al., 1989; Chafee and Goldman-Rakic, 1998; Courtney et al., 1998; Mendoza-Halliday et al., 2014). On the other hand, the sensory recruitment hypothesis posits that working memory storage is mediated by the same areas involved in the initial sensory processing of stimuli and that parietal and prefrontal areas instead serve as a source of top-down biasing signals (Pasternak and Greenlee, 2005; D'Esposito and Postle, 2015). Recently, it has been suggested that working memory contents are distributed across number of cerebral cortical areas including both sensory and prefrontal regions (Ester et al., 2015; Lee and Baker, 2016; Serences, 2016; Christophel et al., 2017). Despite the progress that has been made in characterizing the neural substrates of working memory encoding and storage, the involvement of subcortical structures has received little consideration. A number of studies have reported that portions of the cerebellum are actively involved in attention and working memory task performance (Allen et al., 1997; Desmond et al., 1997; Chen and Desmond, 2005a; Kirschen et al., 2005; Hayter et al., 2007; Baier et al., 2010; Stoodley et al., 2012; Striemer et al., 2015; Brissenden et al., 2016, 2018; Guell et al., 2018a; King et al., 2019). However, as these studies either did not include a VWM maintenance component or employed a block design in which sample, delay, and probe responses could not be distinguished, it is unclear whether the cerebellar contribution to these tasks reflects working memory storage processes or rather a more generalized attentional control mechanism. The cerebellum has been proposed to play a predictive role in working memory (Sheu et al., 2019), mirroring cerebellar models of motor control (Wolpert et al., 1998). In order to support a predictive function in WM, the cerebellum would need to represent the item-specific contents of working memory. In the current study, we sought to directly investigate whether specific portions of the cerebellum contribute to the encoding and maintenance of mnemonic representations.
Early research aimed at determining where working memory contents are stored in the brain was primarily informed by either lesion-induced deficits in working memory task performance or the observation of sustained, elevated activity over extended delays (Funahashi et al., 1989; Courtney et al., 1998; Postle et al., 2000; Wager and Smith, 2003; Sreenivasan and D'Esposito, 2019). However, working memory deficits following lesioning or elevated delay period activity could simply indicate that an area influences storage in other areas and does not imply that an area contains a representation of working memory contents (D'Esposito and Postle, 2015; Lee and Baker, 2016; Serences, 2016; Christophel et al., 2017). More recently, researchers have employed multivariate methods to determine whether an area exhibits activity patterns that are selective for the specific stimulus stored in working memory (Serences et al., 2009; Harrison and Tong, 2009; Christophel et al., 2012, 2015, 2018; Riggall and Postle, 2012; Lee et al., 2013; Pratte and Tong, 2014; Ester et al., 2015). As a result, content specificity is now argued to be the critical marker of information storage in VWM (Christophel et al., 2017). Here, we used an event-related delayed recall paradigm and a multivariate encoding model to investigate whether this distributed network extends to the cerebellum.
Materials and Methods
Participants
A total of 16 healthy adult volunteers (nine female) participated in this study. The Institutional Board of Boston University approved the study. All subjects gave written informed consent and were compensated to participate in this study. Subjects were recruited from Boston University and the Greater Boston area. Our sample included six laboratory members (including one author), each of which had extensive experience as an experimental subject in neuroimaging studies. The remaining participants were not lab members but had prior experience performing psychophysical tasks in an fMRI scanner. All subjects possessed normal or corrected-to-normal vision. All subjects completed a ∼30-min behavioral session before being scanned to familiarize them with the task. One subject's behavioral accuracy was found to not significantly differ from chance performance and was excluded from further analysis, leaving us with 15 subjects (nine female).
Magnetic resonance image acquisition
Data were acquired from a 3 Tesla Siemens Prisma scanner located at the Cognitive Neuroimaging Center at Boston University in Boston, Massachusetts using a 64-channel head coil. High-resolution T1-weighted multiecho MPRAGE (sagittal; TR = 2780 ms; TE = 1.32, 3.19, 5.11, 7.03 ms; FA = 7°; 0.8-mm isotropic voxels; 224 slices; FOV = 256 × 256 × 180 mm; in-plane GRAPPA acceleration 2) and T2-weighted (sagittal; TR = 3200 ms; TE = 564 ms; 0.8-mm isotropic voxels; 224 slices; FOV = 256 × 256 × 180 mm; in-plane GRAPPA acceleration 2) structural images were acquired. Functional data were acquired using a multiband gradient-echo echoplanar pulse sequence (Feinberg et al., 2010; Moeller et al., 2010; Xu et al., 2013) with the following acquisition parameters: TR = 2000 ms; TE = 35 ms; FA = 80°; 2.2-mm isotropic voxels; FOV = 207 × 207 × 152 mm; multiband factor (SMS) = 3. A total of 69 slices were acquired with 0% skip, fully covering the cerebral cortex and cerebellum. Image reconstruction was performed using the split-slice GRAPPA method, which has been shown to dramatically reduce the incidence of false positive activation because of signal leakage (Cauley et al., 2014; Todd et al., 2016). Spin echo field maps were also acquired with opposite phase encoding directions (anterior-to-posterior; posterior-to-anterior) and matching parameters to the gradient-echo EPI fMRI timeseries for subsequent EPI distortion correction.
Visual stimuli and experimental paradigm
Stimuli were generated and presented using Python with the PsychoPy software package (Peirce, 2007, 2008) and were projected with a PROPixx digital light processing LED projector (VPixx Technologies Inc.) onto a screen within the scanner bore. Each trial consisted of a trial start cue period (1 s), sample period (1 s), retro-cue and delay period (10 s), and probe period (5 s; Fig. 1A). The beginning of each trial was denoted by the white fixation cross briefly changing to red for 500 ms. A 500-ms fixation interval separated the trial start cue and the sample period. During the sample period, subjects were presented with two circular random dot kinematograms (one in each hemifield) each spanning 8° in diameter. Dots within each patch were 0.05° in size and possessed a density of ∼2 dots per square degree. The dot motion patches were centered 6° to the left or right from fixation along the horizontal meridian. Dots moved with 100% coherence at 1.5°/s and had a limited lifetime of 10 frames (167 ms). The direction of coherent motion on a trial within each patch was drawn from a uniform distribution over 0−315° in 45° increments. A small ±1−10° angular jitter was added randomly to the motion direction on each trial to minimize verbal coding of stimuli; 250 ms following the offset of the sample, an arrow appeared pointing either to the left or right. This “retro-cue” indicated which of the two motion directions subjects should maintain over the subsequent 9.75-s delay period. The use of a retro-cue was a key element of the experimental design as it permitted us to distinguish between initial stimulus encoding and maintenance in VWM by comparing delay period activity between the remembered and non-remembered stimuli. On 20% of trials, an “X” was presented indicating that both items should be dropped. Following the delay period, subjects were presented with a probe display (5 s). The probe display consisted of a circle (8° diameter) centered at fixation and a line extending from the center of the display to a location on the circumference of the circle. Subjects were instructed to rotate the line so that it aligned with the remembered direction of motion. Subjects pressed one key to incrementally rotate the line 10° clockwise, another key to rotate the line 10° counterclockwise, and a third key to flip the line 180°. The initial direction of the line segment was selected randomly among the 36 possible offsets that are an integer multiple of 10° of the to-be-remembered motion direction. The randomization of the probe stimulus ensured that subjects were unable to prospectively plan a specific motor action (i.e., direction of rotation) during the delay period. On drop trials, subjects were instructed to make a random number of button-press responses. Trials were separated by a 7-s intertrial interval. Each fMRI run comprised 20 trials (16 store trials + 4 drop trials; 482 s total). Subjects completed eight runs during the session. In-scanner eye position tracking was planned, but because of hardware issues was not possible during this experiment. Stimulus direction and location (left or right visual hemifield) of the cued stimulus were fully crossed within a run. Additionally, the motion directions of cued and non-cued stimuli were fully crossed across all eight runs.
fMRI data preprocessing and regions of interest (ROIs)
Functional task data first underwent the Human Connectome Project (HCP)'s “minimal” preprocessing pipeline, which comprises gradient nonlinearity distortion correction, motion correction, EPI image distortion correction, co-registration with the subject's T1-weighted image, and nonlinear transformation to MNI152 space (Glasser et al., 2013). The transforms involved in each step of this pipeline were concatenated into a single nonlinear transformation and performed as a single resampling step to reduce interpolation related blurring (Glasser et al., 2013). Voxel timeseries within each run were then standardized and detrended with a third order polynomial. Multivariate analyses were performed on the average response of individual voxels across three TRs beginning 6, 8, and 10 s (TRs 4, 5, and 6) after the start of each trial, consistent with previous working memory decoding studies (Harrison and Tong, 2009; Ester et al., 2015). Additionally, we generated stimulus reconstructions across time by applying our decoding model to a sliding window of the average response across three consecutive TRs (0–4, 2–6 s, etc.)
Resting-state data underwent additional preprocessing using custom scripts in MATLAB. The following preprocessing steps were performed: linear interpolation across high-motion time-points (>0.5 mm FD; Power et al., 2012; Carp, 2013), application of a fourth-order Butterworth temporal bandpass filter to extract frequencies between 0.009 and 0.08 Hz, mean “grayordinate” signal regression (MGSR; Burgess et al., 2016) and censoring of high-motion time-points (Power et al., 2012).
ROIs were defined from publicly available atlases. Cortical ROIs corresponding with a combined IPS0-IPS3 ROI [left (L): 670 vertices; right (R): 606 vertices] and a combined V1-V3 ROI (L: 1008 vertices; R: 960 vertices) were defined in each hemisphere from a probabilistic retinotopy atlas (Wang et al., 2015). A putative VWM region in dorsolateral frontal cortex was defined from the 6a label from the Glasser multimodal parcellation in both hemispheres (L: 360 vertices; R: 339 vertices; Glasser et al., 2016); area 6a corresponds with the intersection of the superior frontal sulcus and the superior branch of the precentral sulcus, and exhibits robust working memory activation in the HCP VWM (2-back vs 0-back) task contrast (Glasser et al., 2016). This ROI will be referred to as superior frontal junction (SFJ) throughout the remainder of the manuscript. The cerebellar spatially unbiased infratentorial template (SUIT) atlas was used to define an ROI for each lobule within the cerebellum in each hemisphere (L: 133–1866 voxels; R: 140–1765 voxels; Diedrichsen et al., 2009). To exclude any potential signal spill-over (resulting from the nonlinear spatial transformation procedure) from ventral visual cortex into cerebellum (Buckner et al., 2011; van Es et al., 2019), we applied a dilated (6 mm) cortical ribbon mask to dorsal cerebellar lobule ROIs (lobules I–IV, V, and VI and Crus I). As our prior work indicates that an area spanning both lobule VIIb and lobule VIIIa is recruited by attention and working memory tasks along with fronto-parietal areas (Brissenden et al., 2016, 2018), we created a combined lobule VIIb/VIIIa ROI (L: 1660 voxels; R: 1731 voxels). An oculomotor vermis ROI (L: 99 voxels; R: 55 voxels) was defined using a publicly available cerebellar retinotopy atlas produced by an analysis of the HCP retinotopy dataset (van Es et al., 2019). We further created ROIs subdividing VIIb/VIIIa using resting-state functional connectivity. To do so, a resting-state correlation with each cortical ROI (V1-3, IPS0-3, SFJ) was computed for each voxel within lobule VIIb/VIIIa. As each cerebral cortical hemisphere is anatomically connected with the contralateral hemisphere of cerebellar cortex (Strick, 1985; Schmahmann and Pandya, 1997), left hemisphere cortical seeds were correlated with right hemisphere lobule VIIb/VIIIa voxels and vice versa. We then created ROIs for each subject corresponding with the top 25% of correlations within lobule VIIb/VIIIa for each cortical seed. A percentile threshold was used to account for potential SNR differences between subjects and to ensure that ROIs were of similar size across subjects. This analysis revealed a differential pattern of connectivity within lobule VIIb/VIIIa for each cortical seed. To quantify the location of each cortical area's connectivity within lobule VIIb/VIIIa across subjects, we summed 25th percentile ROIs across subjects to create a probabilistic subject count map for each cortical ROI. The top 25% of resting-state correlations for IPS0-3 and SFJ were found to be concentrated in the medial and lateral portions of lobule VIIb/VIIIa for the majority of subjects. In contrast, the top 25% of V1-3 correlations was more evenly distributed throughout lobule VIIb/VIIIa across subjects. We then subdivided lobule VIIb/VIIIa on the basis of this differential connectivity. Medial (L: 376 voxels; R: 322 voxels) and lateral (L: 105 voxels; R: 128 voxels) VIIb/VIIIa ROIs were defined by taking the maximum across the IPS0-3 and SFJ subject count maps and then thresholding at six or more (i.e., >33.3% of) subjects. A central lobule VIIb/VIIIa ROI (L: 161 voxels; R: 140 voxels) intermediate to the medial and lateral ROIs was defined by thresholding the V1-3 subject count map at six or more subjects and then excluding any voxel that already belonged to either the medial or lateral ROI.
Statistical analysis
Behavioral analysis
To assess each participant's behavioral accuracy, we computed the mean absolute value of the angular difference between reported and actual motion directions. We further fit several different mixture models to each participant's response error distribution using the MemToolBox package (Suchow et al., 2013). We fit the following models: standard mixture model (Zhang and Luck, 2008), swap model (Bayset al., 2009), variable precision model (van den Berg et al., 2012; Fougnie et al., 2012), and a swap + variable precision model (van den Berg et al., 2014). The standard mixture model fits participants' response errors with a mixture of two distributions, a uniform distribution and a von Mises distribution, reflecting (1) trials in which items were not stored in working memory and (2) the precision of items successfully represented in working memory [parameters: guess rate (g) and precision (SD)]. The swap model is an extension of the standard mixture model that includes an additional parameter representing the proportion of trials for which participants incorrectly reported the identity of the non-cued item [parameters: g, SD, swap rate (b)]. The variable precision model, on the other hand, assumes that the precision of a participant's reports is itself drawn from a higher-order variability distribution [parameters: g, mean precision (SDµ), precision variability (SDσ)]. The swap + variable precision combines swap and variable precision models and includes parameters associated with both non-target reports and precision variability (parameters: g, b, SDµ, and SDσ). Each model additionally included a bias parameter to allow the central tendency of the fitted distribution to differ from zero. To assess the relative goodness of fit between models, we computed the Akaike Information Criterion (AIC; Akaike, 1974) for each model:
where L is the model log-likelihood and k is the number of free parameters. As AIC assumes infinite data and will tend to overfit for smaller datasets, we computed a corrected version of AIC (AICc; Hurvich and Tsai, 1989):
where n represents the number of working memory trials used to fit the model. The variable precision model was determined to best fit our data (mean AICc difference ± SEM: 20.17 ± 5.29 vs standard mixture model, 22.25 ± 5.27 vs swap model, 2.43 ± 0.47 vs swap + variable precision model). Consequently, we report mean parameters for the best fitting variable precision model across subjects associated with guess rate (g), mean precision (SDµ), and precision variability (SDσ).
Generative model
We employed a generative model-based approach developed by van Bergen et al. (2015). This model assumes voxel responses can be modeled as the linear weighted sum of 8 motion direction selective neuronal populations or channels. Each channel was represented as a half-wave rectified (co)sinusoid raised to the fifth power:
where θk indicates the preferred direction of the kth population or channel. Channel basis functions were maximally tuned at one of eight equally spaced angles (0°, 45°, 90°, 135°, 180°, 225°, 270°, or 315°).
Model fitting and assessment was performed using an eightfold leave-one-run-out cross-validation scheme, which iteratively partitioned the data into a training set (B1) and a test set (B2). The training set was expressed as an n × m matrix B1, where n is the number of trials in the training set and m is the number of voxels. To account for retinotopic preferences, the training and test sets were restricted to trials in which the stimulus appeared in the ROI's preferred hemifield (contralateral for cortical ROIs and ipsilateral for cerebellar ROIs; Brissenden et al., 2018; van Es et al., 2019). An n × k hypothetical channel output matrix C1 (where k is the number of channels) represented the idealized channel responses for each trial given the remembered motion direction. A k × m channel weight matrix W related the observed voxel responses to the hypothetical channel responses via a general linear model of the form:
The weight matrix was estimated via ordinary least squares as follows:
Then for each trial in the test set, we used Bayes rule to obtain a posterior probability distribution over motion direction indicating which stimulus value was most probable given the observed pattern of BOLD responses.
where s denotes a particular stimulus value and b2 indicates a single trial/row of test set BOLD matrix (B2). The prior, p(s), was set to 1 for all angles as each stimulus value was presented equally often and we did not want to bias decoding toward particular motion directions (van Bergen et al., 2015). To avoid numerical underflow, we computed log-probability rather than probability. The conditional log-probability of a voxel activation pattern given a specific stimulus value was defined as:
where c2 denotes a vector of hypothetical channel responses corresponding with the stimulus value s and Σ denotes the m × m voxel covariance matrix estimated from the training set. As computing the above conditional probability requires computing the inverse of the voxel covariance matrix, we performed a procedure described in Naselaris et al. (2009) if Σ was singular. First, we performed principal component analysis on the predicted responses (
) for the training set. p components were chosen with the constraint that each component explained at least 5% of the total variance. Using the resulting m × p projection matrix, we then projected the predicted (
) and actual (B1) responses from the training set onto the first p principal components. These dimensionality-reduced responses were then normalized to unit length:
where P represents the computed projection matrix, 
and 
respectively denote a single row from the predicted (
) and actual (B1) response matrix, and 
and 
, respectively, denote a single row from the dimensionality-reduced n × p predicted (
) and actual (
response matrix. The voxel covariance matrix was then computed as:
where Σ* is a p × p matrix. During the testing phase, the same procedure (projection onto the first p components using the projection matrix obtained from the training set followed by normalization) was performed on the predicted and actual responses for the test set. Thus, the conditional log-probability of the dimensionality-reduced pattern of voxel responses given a specific stimulus value was re-expressed as:
To assess whether an area contains stimulus-specific information across trials, we then circularly shifted trial log-posterior distributions so that the center value corresponded with the presented stimulus on that trial. Probability densities, as normalized ratios, must be averaged using the geometric mean rather than the arithmetic mean (Fleming and Wallace, 1986; Nelson, 2017). This is equivalent to computing the arithmetic mean of log-probability and then exponentiating the result. Consequently, we averaged the shifted log-probability distributions across trials to produce a subject-level log-posterior probability distribution. We then averaged subject log-posterior distributions to form a group-level log-posterior probability distribution. As the training and testing procedure was performed separately for each hemisphere, a group-level log-posterior probability distribution was first generated for each hemisphere and then these distributions were averaged across hemispheres. The bilateral group-level log-posterior probability distribution was then exponentiated and divided by a normalization constant to compute the geometric mean and yield a probability distribution. To assess the degree to which ROI activity patterns contained information about the presented stimulus, we performed non-parametric permutation tests comparing the probability of the presented stimulus [corresponding with p(0°|b) for circularly shifted distributions] to an empirical null distribution (Golland and Fischl, 2003). To estimate a null distribution, we permuted the motion direction labels within each run before model fitting and testing. This procedure was repeated 1000 times for each subject, yielding a group-level probability distribution for each permutation. We then compared the empirical null distribution to the probability obtained with the true class labeling (α < 0.05, one-sided). A p value was computed as [(# of permutation probabilities ≥ actual probability) + 1]/(N permutations + 1) (Phipson and Smyth, 2010). Correction for multiple comparisons across ROIs was performed using the Holm–Bonferroni procedure (Holm, 1979). Inspired by cluster correction approaches for group-level fMRI analyses (Heller et al., 2006; Woo et al., 2014; Eklund et al., 2016), time-series decoding comparisons were corrected using a nonparametric cluster correction procedure to account for the dependency across sliding TR windows. A cluster-defining threshold was set at p < 0.05. Adjacent timepoints with below-threshold p values (<0.05) were considered a cluster. Permutation probabilities were also converted to p values by comparing each permutation to the remaining permutations and the actual probability, i.e., count number of times permuted probabilities or the actual probability is greater than or equal to the current permutation and then divide by the number of permutations plus 1. We then created a maximum null cluster size distribution by thresholding and recording the maximum cluster size for each permutation across TR windows (Nichols and Holmes, 2002). A p value for the observed cluster size was computed as [(# of permutation cluster sizes ≥ actual cluster size) + 1]/(N permutations + 1) (Phipson and Smyth, 2010). It should be noted that permutation tests impose a minimum p value that is determined by the number of permutations performed. Here, the minimum possible p value before correction for multiple comparisons was 1/1001 = 0.000999. The minimum p value was greater following multiple comparison correction. We additionally examined whether the posterior probability of the presented stimulus given an area's activity patterns was correlated between ROIs across trials. To do so, we computed pairwise correlations between trial-wise estimates of the posterior probability of the presented stimulus for lobule VIIb/VIIIa, IPS0-3, SFJ, and V1-3. This analysis can be likened to a β-series correlation analysis (Rissman et al., 2004) that examines the correspondence between areas in terms of representational content rather than task-evoked BOLD response. Permutation tests were performed by repeatedly shuffling trial-wise posterior probabilities within each subject before computing the correlation between areas. We further investigated the relationship between area decoding correlations for each subject and parameters derived from the variable precision mixture model fitted to each subject's response error distribution. We again performed permutation tests by independently shuffling subject-wise decoding correlations and model parameters before computing the correlation between the two.
Results
To investigate VWM stimulus specificity in the cerebellum, we had participants perform a delayed motion direction recall task while measuring BOLD responses using fMRI (Fig. 1A). On each trial, participants were simultaneously presented with two patches of dot motion (100% coherence), presented to the right and left of fixation. After stimulus offset, participants were retroactively cued to maintain the direction of motion of one of the patches over a blank delay interval (9.75 s) or drop both items. Following the delay period, participants then adjusted a probe stimulus to match the remembered direction of motion. Behavioral performance was assessed by computing a distribution of recall errors over trials for each subject (Fig. 1B). Trial errors were found to be clustered around the cued direction of motion, indicating that participants were able to successfully maintain the direction of motion of the cued item. The average mean absolute value of the angular difference between reported and actual motion directions was 16.75 ± 1.09° SEM. We additionally fit a computational model to participants' response error distribution that enabled us to estimate for each subject the proportion of trials in which they guessed, the mean precision with which they remembered stimuli, and the degree to which precision was variable across trials (see methods; van den Berg et al., 2012; Fougnie, Suchow and Alvarez, 2012). The average guess rate (g) across subjects was 9.29 ± 2.64% SEM. The average precision (SDµ) across subjects was 8.04 ± 1.41° and the average standard deviation of precision (SDσ) was 16.67 ± 2.28°.
Behavioral task paradigm and behavioral performance. A, Participants were presented with two patches of moving dots. 250 ms following the offset of the stimulus presentation period participants were postcued to remember the direction of motion of the moving dots. Following a 10-s delay, participants were given 5 s to adjust an oriented line segment to match the remembered direction of motion with key presses. The initial orientation of the line segment was randomized with respect to the cued direction of motion on each trial. B, Kernel density estimates of recall error for all participants. The colored lines denote density estimates of recall error for each individual participant and the black line reflects the density estimate of recall error across all participants.
To ensure that participants were sensitive to the actual stimulus direction on each trial and not simply storing a categorical representation (i.e., 0°, 45°, 90°, 135°, 180°, 225°, 270°, or 315°), we investigated the relationship between the directional jitter applied on each trial (±1–10°) and response error. To do so, we examined recall errors for five categories of trials based on jitter direction and size [zero jitter, small clockwise jitter, small counterclockwise jitter (<5°), large clockwise jitter, and large counterclockwise jitter (>5°)]. We did not find an overall effect of jitter size/direction on signed error, F(4,1912) = 1.96, p = 0.097. Pairwise comparisons between each jitter size/direction and zero jitter further found no significant differences (all p > 0.37), indicating that participants did not exhibit any directional bias in their responses relative to non-jittered trials. We additionally found no differences in absolute error between jitter conditions (all p > 0.43), which showed that recall precision did not depend on the direction or size of jitter applied and therefore that subjects were encoding the actual direction of motion rather than a categorical representation of the non-jittered direction.
Stimulus encoding model
To quantify feature selectivity during working memory, we employed a generative model-based decoding approach. This analysis uses information about voxel motion direction preferences to produce a posterior probability distribution over motion direction for each trial given the observed pattern of voxel responses during the delay period (see methods). To produce a group-level reconstruction of the probability of a particular stimulus given the observed responses, we circularly shifted trial posterior probability distributions so 0° corresponded with the presented motion direction and then averaged the shifted distributions across trials and participants. If an ROI contains feature-selective information about the remembered direction of motion then this averaged posterior distribution should peak at or near 0°. We first present our decoding results for the last three TRs of the delay period (TRs beginning 6–10 s after trial onset), consistent with previous working memory decoding studies (Harrison and Tong, 2009; Ester et al., 2015). Centered and averaged posterior distributions for V1-3, IPS0-3, and SFJ are shown in Figure 2. Each cortical ROI was found to contain a robust representation of the direction of motion maintained in working memory (IPS0-3: p = 0.004 corrected; SFJ: p = 0.004 corrected; V1-3: p = 0.004 corrected; minimum possible p value following correction was 0.004; Fig. 2A–C). The representation of the non-remembered stimulus in these areas also appears to persist into the delay period but in a weaker form (IPS0-3: p = 0.004 corrected, SFJ: p = 0.044 corrected; V1-3: p = 0.027 corrected; Fig. 2D–F). The difference between remembered and non-remembered stimuli for the average across the last three TRs of the delay period is significant for SFJ (p = 0.027 corrected). Although numerically greater the difference for IPS0-3 and V1-3 does not reach significance (p = 0.131 uncorrected and p = 0.129 uncorrected). However, a difference does emerge at the following TR window (average across TRs 5, 6, and 7; IPS: p = 0.079 corrected; V1-3: p = 0.002 corrected). These results replicate prior reports of working memory content specificity in visual, parietal and frontal cortex (Christophel et al., 2012, 2017; Riggall and Postle, 2012; Lee et al., 2013, 2013; Sprague et al., 2014; Ester et al., 2015).
Group-level cortical encoding model results. Shifted and averaged posterior probability distributions over motion direction indicating which stimulus is most probable given the observed pattern of vertex responses during the delay period for (A) IPS0-3, (B) SFJ, and (C) V1-3. Posterior probability over motion direction for non-remembered stimuli for (D) IPS0-3, (E) SFJ, and (F) V1-3. Shaded ribbon reflects a bootstrap estimate of SEM.
We next examined whether cerebellar lobule VIIb/VIIIa delay-period responses were also selective for a stimulus stored in working memory. lobule VIIb/VIIIa was preselected for analysis as our prior work reported VWM task activation localized to this portion of cerebellum (Brissenden et al., 2016, 2018). A group-level reconstruction of the probability of a stimulus given the pattern of voxel responses in lobule VIIb/VIIIa was found to peak at the remembered direction of motion (p = 0.004 corrected; Fig. 3A). We also found a significant, albeit weak, representation of the non-remembered stimulus in lobule VIIb/VIIIa (p = 0.039 uncorrected; Fig. 3B). To test whether the strength of representation was stronger for remembered or non-remembered stimuli we compared the observed difference in probability for remembered and non-remembered items to a null probability difference distribution. The probability of the remembered stimulus given the observed pattern of lobule VIIb/VIIIa voxel responses was significantly greater than the probability of the non-remembered stimulus (p = 0.0025 corrected). This difference indicates that the representation of the remembered stimulus in VIIb/VIIIa during the delay period is working memory specific and cannot be explained by residual stimulus-driven responses. We also examined the hemispheric specificity of working memory representations in lobule VIIb/VIIIa by performing the same analysis for stimuli presented in the non-preferred hemifield (contralateral hemifield for cerebellum). We found a significant representation of the remembered item in the contralateral hemifield in lobule VIIb/VIIIa (p = 0.004). However, this representation did not significantly differ from the non-cued item (p = 0.788). Thus, in contrast to items presented in the preferred hemifield (ipsilateral) we are unable to claim that the VIIb/VIIIa representation of items in the non-preferred hemifield is specific to VWM.
Lobule VIIb/VIIIa encoding model results. A, Posterior probability over motion direction for remembered stimuli. B, Posterior probability over motion direction for non-remembered stimuli. Shaded ribbon reflects bootstrap SEM. C, Time-resolved encoding model results. Each row of the heat map represents an averaged and shifted probability distribution generated by applying an encoding model to an average of three consecutive TRs. tr123 – TRs 1, 2, 3; tr234 – TRs 2, 3, 4; tr345 – TRs 3, 4, 5; tr456 – TRs 4, 5, 6; tr567 – TRs 5, 6, 7; tr678 – TRs 6, 7, 8.
We next investigated the relationship between decoding in lobule VIIb/VIIIa and cerebral cortical regions. To do so, we extracted the probability of the presented stimulus from shifted posterior distributions for each trial [i.e., p(0°|b) for circularly shifted distributions] and then correlated these values between areas for each subject. lobule VIIb/VIIIa decoding was positively correlated with decoding in each cortical ROI and permutation tests indicated that correlations averaged across subjects were greater than what would be expected by chance (VIIb/VIIIa–IPS0-3: r = 0.171 ± 0.052, p = 0.006 corrected; VIIb/VIIIa–SFJ: r = 0.142 ± 0.050, p = 0.006 corrected; VIIb/VIIIa–V1-3: r = 0.148 ± 0.038, p = 0.006 corrected; 0.006 was the minimum possible corrected p value for 6 comparisons). Furthermore, the observed trial-wise correlations between lobule VIIb/VIIIa and cortical ROIs were comparable to those between cortical ROIs (IPS0-3–SFJ: r = 0.281 ± 0.037, p = 0.006 corrected; IPS0-3–V1-3: r = 0.161 ± 0.039, p = 0.006 corrected; SFJ–V1-3: r = 0.101 ± 0.025, p = 0.006 corrected). Taken together, these results indicate that fluctuations in information content within lobule VIIb/VIIIa and canonical VWM areas are linked on trial-wise basis.
We also found that the correlation in decoding between certain areas predicted behavioral recall precision. The correlation between IPS0-3 and lobule VIIb/VIIIa and the correlation between IPS0-3 and SFJ for each participant were found to be significantly correlated with the mean precision parameter from the variable precision model fitted to each participant's error distribution [cor(IPS0-3–VIIb/VIIIa, SDµ): r = −0.59, p = 0.042 corrected; cor(IPS0-3–SFJ, SDµ): r = −0.61, p = 0.04995 corrected]. No significant correlations were identified between the other model parameters (g and SDσ) and the correlation in decoding between areas. These results indicate that subjects for which trial-by-trial fluctuations in information content were more tightly linked between IPS and both cerebellum and frontal cortex remembered stimuli with greater precision.
Time-resolved lobule VIIb/VIIIa stimulus probability reconstructions for both remembered and non-remembered items were produced by applying the encoding model to a sliding window consisting of three consecutive TRs (Fig. 3C). A non-parametric cluster-based correction procedure, a method commonly used to correct for multiple correlated comparisons in group-level fMRI analyses (Heller et al., 2006; Woo et al., 2014; Eklund et al., 2016), revealed a significant effect from the second to fifth TR window for remembered items (cluster p = 0.003). On the other hand, the non-remembered item was not significantly represented across time (p = 0.224). Furthermore, the difference in probability between remembered and non-remembered items was significantly different from the second to fifth TR window (p = 0.002). Thus, cerebellar lobule VIIb/VIIIa delay period responses appear to be content-specific and this specificity is sustained throughout the delay period.
Working memory stimulus specificity was further shown to be restricted to lobule VIIb/VIIIa within the cerebellum. We were unable to reconstruct the remembered stimulus from any other cerebellar lobule (all p > 0.05 uncorrected), with the exception of Crus II (p = 0.027 uncorrected), which has been previously implicated in working memory (Guell et al., 2018a), exhibits robust functional connectivity with the cognitive control network (Buckner et al., 2011; Guell et al., 2018b), and is adjacent to lobule VIIb/VIIIa (Fig. 4). However, Crus II only weakly represented the remembered stimulus and this effect did not survive correction for multiple comparisons (p = 0.216 corrected). We additionally examined whether an oculomotor vermis ROI, an area known to be highly responsive to saccades (Noda and Fujikado, 1987; Ohtsuka and Noda, 1995; Takagi et al., 1998; Thier et al., 2002; Soetedjo and Fuchs, 2006), could decode the remembered stimulus to determine whether saccadic eye movements could explain the observed stimulus selective activity. We were unable to recover a significant representation of the presented stimulus from oculomotor vermis delay period activity (p = 0.098 uncorrected; Fig. 5). Collectively, these results indicate that the encoding of stimulus-specific mnemonic representations is limited to lobule VIIb/VIIIa within the cerebellum. Furthermore, the lack of evidence for feature-selective representations elsewhere in the cerebellum argues against a purely motor explanation for lobule VIIb/VIIIa content specificity, as motor processes would be expected to drive areas of the cerebellum previously implicated in motor functions such as the anterior lobe (lobules I–V) and the oculomotor vermis.
Cerebellar lobule encoding model results. Group-level posterior probability distributions for cerebellar lobules outside of Lobule VIIb/VIIIa. Inset shows flatmap representation of cerebellum (Diedrichsen and Zotow, 2015). Colors and dotted lines demarcate cerebellar lobules. Shaded ribbon reflects bootstrap SEM.
Oculomotor vermis encoding model results. Averaged and shifted group-level posterior probability distribution for oculomotor vermis. Inset depicts oculomotor vermis ROI overlaid on cerebellar flatmap. Shaded ribbon reflects bootstrap SEM.
We next sought to determine whether a particular subdivision of lobule VIIb/VIIIa was driving the observed working memory stimulus specificity. To do so, we performed a seed-based resting-state functional connectivity analysis using V1-3, IPS0-3, and SFJ ROIs as seeds. As each cerebral cortical hemisphere is preferentially connected with the contralateral hemisphere of cerebellum, these analyses were performed separately for each hemisphere (right hemisphere cerebral cortical ROI to left hemisphere cerebellum and vice versa). For each subject, we then identified the top 25% of correlations within lobule VIIb/VIIIa for each cortical ROI. We then computed the overlap across subjects by summing individual percentile ROIs across subjects for each resting-state analysis (Fig. 6). Cortical ROIs exhibited differential connectivity across lobule VIIb/VIIIa. IPS0-3 was preferentially connected with the medial portion of lobule VIIb across subjects, and to a lesser extent the lateral portion of lobule VIIb/VIIIa (Fig. 6A), while largely sparing the central portion. SFJ connectivity, on the other hand, was equally concentrated in medial VIIb and lateral VIIb/VIIIa, again sparing the central portion, across subjects (Fig. 6B). V1-3 connectivity was more evenly distributed across all of lobule VIIb/VIIIa across subjects (Fig. 6C). Using these subject count maps, we then subdivided lobule VIIb/VIIIa to create ROIs localized to the medial, central, and lateral portions of lobule VIIb/VIIIa bilaterally (see Materials and Methods; Fig. 7A). We then applied the encoding model to each of these ROIs' delay-period activity to determine which portion of lobule VIIb/VIIIa exhibited feature-selection activation. Medial lobule VIIb/VIIIa was shown to encode feature-selective information about the remembered stimulus (p = 0.004 corrected; Fig. 7B), but little information about the non-remembered stimulus (p = 0.09 uncorrected; Fig. 7C). Additionally, medial lobule VIIb/VIIIa exhibited a significant difference in decoding between remembered and non-remembered stimuli (p = 0.02 corrected). Lateral lobule VIIb/VIIIa also exhibited significant decoding of the remembered stimulus (p = 0.004 corrected; Fig. 7F), but was also shown to significantly decode the non-remembered stimulus (p = 0.036 corrected; Fig. 7G). The difference in decoded probability between remembered and non-remembered stimuli was not significant (p = 0.169 uncorrected), indicating that the stimulus representation in lateral lobule VIIb/VIIIa may be in part driven by a perceptual trace rather than a working memory trace. Nevertheless, combining the medial and lateral VIIb/VIIIa ROIs into a single ROI yielded significant decoding of the remembered stimulus (p = 0.004 corrected; Fig. 7H), but not the non-remembered stimulus (p = 0.122 uncorrected; Fig. 7I), and an even more robust difference between remembered and non-remembered stimuli (p = 0.0002 corrected). In contrast, the central portion in each hemisphere of lobule VIIb/VIIIa did not exhibit evidence for a representation of either the remembered or non-remembered stimulus (remembered: p = 0.882 uncorrected; non-remembered: p = 0.111 uncorrected; difference: p = 0.953 uncorrected; Fig. 7D,E). In sum, portions of lobule VIIb/VIIIa that are preferentially connected with fronto-parietal cortex exhibit content-specific activation during the delay-period, while portions of lobule VIIb/VIIIa that are preferentially connected with visual cortex do not.
Cerebellar Lobule VIIb/VIIIa functional connectivity with cortical ROIs. The top 25% of correlations across Lobule VIIb/VIIIa with either (A) IPS0-3, (B) SFJ, or (C) V1-3 were identified for each participant. Values within Lobule VIIb/VIIIa reflect the number of subjects for which that vertex is included in the top 25% of correlations.
Cerebellar Lobule VIIb/VIIIa subdivision encoding model results. A, Medial, central, and lateral VIIb/VIIIa ROIs defined using probabilistic functional connectivity subject count maps depicted in Figure 6. B, C, Medial Lobule VIIb/VIIIa ROI group-level posterior probability distribution for (B) remembered and (C) non-remembered items. D, E, Central Lobule VIIb/VIIIa posterior distribution for (D) remembered and (E) non-remembered items. F, G, Lateral Lobule VIIb/VIIIa posterior distribution for (F) remembered and (G) non-remembered items. H, I, Combined medial and lateral Lobule VIIb/VIIIa posterior distribution for (H) remembered and (I) non-remembered items. Shaded ribbon reflects bootstrap SEM.
Discussion
The findings presented here demonstrate for the first time that cerebellar lobule VIIb/VIIIa encodes VWM representations. Prior research had shown that lobule VIIb/VIIIa is recruited by VWM and attention tasks (Chen and Desmond, 2005b; Kirschen et al., 2005; Brissenden et al., 2016, 2018), yet the extent to which lobule VIIb/VIIIa is involved in working memory storage processes or rather more generalized attentional processes was previously unclear. Here, we used a retro-cue paradigm to examine working memory storage independent of perception-related, attention-related, and retrieval-related processes. An encoding model of motion direction revealed that lobule VIIb/VIIIa robustly represented the remembered stimulus. Moreover, the representation of the remembered stimulus was significantly stronger than the representation of the non-remembered stimulus, indicating that this representation can be attributed to a working memory trace rather than a perceptual one. Delay-period lobule VIIb/VIIIa stimulus-selectivity was further shown to be robustly correlated with selectivity in cortical regions previously implicated in VWM. Notably, the correlation in decoding between IPS and lobule VIIb/VIIIa was strongly related to behavioral recall precision. Among cerebellar structures, stimulus representation during the delay period was unique to lobule VIIb/VIIIa. We were unable to recover the remembered direction of motion from other portions of the cerebellum. Further analysis demonstrated that specific subdivisions of lobule VIIb/VIIIa exhibit stimulus-specific responses. Medial VIIb and lateral VIIb/VIIIa were shown to be preferentially connected with frontoparietal cortex, while the intermediate area was found to be preferentially connected with early visual cortex. Subsequent analysis demonstrated that both medial and lateral portions of lobule VIIb/VIIIa encoded remembered items, while the intervening central portion did not. It is important to note that these same medial and lateral portions have been shown to contain topographic representations of the visual field (Brissenden et al., 2018; van Es et al., 2019). Taken together, these findings indicate that stimuli stored in VWM are represented by the cerebellar component of the dorsal attention network (Brissenden et al., 2016).
The current results provide further evidence for cerebellar contributions to cognitive function. It has become increasingly clear that cerebellar function is not limited to the motor domain and that substantial functional heterogeneity exists within cerebellar cortex (Stoodley and Schmahmann, 2009; Buckner et al., 2011; Stoodley et al., 2012; Brissenden et al., 2016, 2018; Guell et al., 2018a,b; van Es et al., 2019; King et al., 2019). Along with recent work demonstrating the existence of visual field representations in lobule VIIb/VIIIa (Brissenden et al., 2018; van Es et al., 2019), this study demonstrates that the cerebellum possesses considerable representational specificity for items within the focus of attention or working memory. It remains to be seen whether areas of the cerebellum associated with other cognitive domains exhibit similar specificity in terms of their representational content.
A number of prior studies have provided evidence for the recruitment of the cerebellum by both visual and verbal working memory (Desmond et al., 1997; Chen and Desmond, 2005b; Kirschen et al., 2005; Hayter et al., 2007; Stoodley and Schmahmann, 2009; Stoodley et al., 2012; Brissenden et al., 2016; Guell et al., 2018a; King et al., 2019). These studies revealed activation across a broad swath of cerebellar cortex including lobule VI, Crus I, Crus II, lobule VIIb, and lobule VIIIa in response to a variety of working memory paradigms. This body of research includes our prior work in which we showed that the portion of lobule VIIb/VIIIa that is functionally connected with cortical dorsal attention network areas responds in a load-dependent manner to a VWM change detection task (Brissenden et al., 2016). However, as the methods used in these prior studies did not lend themselves to revealing content-specific activity, it was previously unclear whether the cerebellum could be specifically implicated in the storage of VWM representations. In the present study, we employed an event-related retro-cue paradigm, which enabled the dissociation of activity related to remembered and non-remembered items during the maintenance period. Consequently, we were able to show that specific cerebellar areas encoded the identity of a stimulus during the delay-period and that this selectivity was substantially more robust for remembered items. These findings extend the prior literature on cerebellar contributions to working memory to provide novel evidence that the cerebellum is actively involved in the persistent storage of information in working memory.
A substantial body of work has attempted to identify the locus or loci of a working memory storage buffer. Our results extend recent proposals that working memory contents are distributed across a number of areas (Lee and Baker, 2016; Serences, 2016; Christophel et al., 2017) by showing that working memory representations are additionally encoded in the cerebellum. Central to the distributed working memory network hypothesis is the notion that storage across different areas reflects different levels of representational abstraction (Christophel et al., 2017). Recent work has proposed a predictive role for cerebellum in verbal working memory (Sheu et al., 2019). Motor control studies suggest that cerebellum serves as a “forward-model” that predicts actions (Wolpert et al., 1998). Given the notable homogeneity of cerebellar circuit architecture, similar computations may be conducted by cognitive portions of the cerebellum (Ramnani, 2006). Follow-up work will need to more finely probe the nature of representations in cerebellum to elucidate its unique role in VWM storage processes.
As the cerebellum has been extensively implicated in oculomotor processes (Ron and Robinson, 1973; Fujikado and Noda, 1987; Voogd and Barmack, 2006; Baier et al., 2009), a candidate alternative explanation for our findings is that VIIb/VIIIa activity patterns reflect eye movement processes rather than a VWM representation. This interpretation is unlikely for several reasons. First, we could not decode remembered motion direction from oculomotor vermis, an area that has been extensively implicated in the control of eye movements (Noda and Fujikado, 1987; Ohtsuka and Noda, 1995; Takagi et al., 1998; Thier et al., 2002; Soetedjo and Fuchs, 2006). Second, participants were experienced observers who extensively practiced maintaining fixation before scanning. Although hardware problems prevented us from tracking eye position in the present study, in our prior work with lateralized VWM targets, subjects exhibited no difference in either eye position or movement between “remember-left” and “remember-right” conditions (Brissenden et al., 2018). Third, our findings cannot be explained by the two most obvious saccade strategies. One strategy would be to make eye movements during stimulus presentation along the direction of motion. However, as participants did not know which dot motion patch would be cued, there is no reason that eye movements would be selective for the direction of motion of the subsequently cued item and not the non-cued item. Another strategy would be to make a saccade to the cued hemifield following the presentation of the cue. Yet, as motion direction and hemifield were dissociated from one another, a saccade toward the cued hemifield would not be predictive of stimulus identity. A less obvious strategy would be to refrain from making eye movements during the presentation of coherent motion and then make a saccade parallel to the direction of motion of the cued item following the presentation of the cue and the offset of the motion stimulus. lobule VIIb/VIIIa stimulus specificity could potentially be explained by this scenario. Concurrent eye-tracking during fMRI scanning would allow future studies to explicitly rule out this possibility.
A recent study in mice reported that the persistent representation of information in frontal cortex depends on the cerebellum (Gao et al., 2018). Head-fixed mice were presented with a sample stimulus that cued one of two actions (left or right lick) to be made following a delay period. Neurons in both frontal cortex and the cerebellar fastigial nucleus were shown to exhibit selectivity for the cued action during the delay period. Critically, optogenetic silencing of fastigial neurons abolished selectivity in frontal cortex and resulted in incorrect behavioral choices. Genetic manipulation of different cell types in the cortico-cerebellar circuit has the potential to further advance our understanding of delay activity in the cerebellum (Low et al., 2018). In the current study, we explicitly dissociated motor planning from VWM storage by randomizing the start position of the probe stimulus relative to the remembered stimulus. As a consequence, subjects were unable to prospectively plan their responses. Thus, our results suggest that the cerebellum is involved in the sustained representation of information in working memory even when stored items are not related with specific actions. Stimulation protocols such as TMS could be used to investigate whether inactivation of specific cerebellar areas in humans similarly interferes with cortical representations of remembered stimuli, as well as behavioral performance (Halko et al., 2014; Esterman et al., 2017).
Cortico-centric models of cognition are pervasive in cognitive neuroscience. Our findings highlight the shortcomings of a narrow focus on cerebral cortex in characterizing the neural mechanisms of working memory storage. The more expansive characterization of the neural substrates of working memory storage suggested by the current results could provide new insights into a wide range of goal-directed behaviors that are known to rely on working memory.
Footnotes
This work was supported by National Institutes of Health Grants R21-EY027703 (to D.C.S.), R01-MH116170 and R01-MH111868 (to M.A.H.), and F31-NS103306 (to S.M.T.), as well as the National Science Foundation Grant DGE-1247312 (to J.A.B.). We thank Dr. Stephanie McMains and Dr. Ross Mair for technical advice.
The authors declare no competing financial interests.
- Correspondence should be addressed to David C. Somers at somers{at}bu.edu or James A. Brissenden at brissend{at}umich.edu
