Abstract
To understand how the brain handles mentally represented information flexibly in the absence of sensory stimulation, working memory (WM) studies have been essential. A seminal finding in monkey research is that neurons in the prefrontal cortex (PFC) retain stimulus-specific information when vibrotactile frequencies were memorized. A direct mapping between monkey studies and human research is still controversial. Although oscillatory signatures, in terms of frequency-dependent parametric beta-band modulation, have been observed recently in human EEG studies, the content specificity of these representations in terms of multivariate pattern analysis has not yet been shown. Here, we used fMRI in combination with multivariate classification techniques to determine which brain regions retain information during WM. In a retro-cue delayed-match-to-sample task, human subjects memorized the frequency of vibrotactile stimulation over a 12 s delay phase. Using an assumption-free whole-brain searchlight approach, we tested with support vector regression which brain regions exhibited multivariate parametric WM codes of the maintained frequencies during the WM delay. Interestingly, our analysis revealed an overlap with regions previously identified in monkeys composed of bilateral premotor cortices, supplementary motor area, and the right inferior frontal gyrus as part of the PFC. Therefore, our results establish a link between the WM codes found in monkeys and those in humans and emphasize the importance of the PFC for information maintenance during WM also in humans.
SIGNIFICANCE STATEMENT Working memory (WM) research in monkeys has identified a network of regions, including prefrontal regions, to code stimulus-specific information when vibrotactile frequencies are memorized. Here, we performed an fMRI study during which human subjects had to memorize vibratory frequencies in parallel to previous monkey research. Using an assumption-free, whole-brain searchlight decoding approach, we identified for the first time regions in the human brain that exhibit multivariate patterns of activity to code the vibratory frequency parametrically during WM. Our results parallel previous monkey findings and show that the supplementary motor area, premotor, and the right prefrontal cortex are involved in vibrotactile WM coding in humans.
Introduction
Flexible handling of mentally represented information is the basis for any higher cognitive process. To investigate how the brain represents mental material in the absence of sensory stimulation, working memory (WM) studies have been essential. For a long time, the dominant view was that major aspects of WM content are maintained in prefrontal regions (Goldman-Rakic, 1995). This view was supported by electrophysiological recordings in monkeys finding that neurons in the prefrontal cortex (PFC) exhibited stimulus-specific firing (Romo et al., 1999). Studies on sensory WM by Romo and colleagues, who explored the retention of vibration frequency of tactile stimuli in nonhuman primates, were particularly influential (Hernández et al., 2000; Romo and Salinas, 2001, 2003a, Romo et al., 2002, 2004). These seminal studies have established direct links between the observed neuronal signals and the actual content that was memorized (Pasternak and Greenlee, 2005), namely, that firing of neurons in the PFC was modulated parametrically depending on the maintained vibration frequency (Romo et al., 1999). Consecutive studies depicted a network of regions exhibiting similar WM codes while controlling for potentially conflating perceptual or decisional processes (Romo and de Lafuente, 2013; Romo et al., 2012). Romo and colleagues were able to show parametric mnemonic representations in secondary somatosensory, medial frontal, and premotor cortices and the PFC (Romo and de Lafuente, 2013; Romo et al., 2012). Likewise, several human oscillatory EEG studies have emphasized the role of the right PFC in maintaining information about the stimulus feature (i.e., vibrotactile frequency; Spitzer and Blankenburg, 2012; Spitzer et al., 2014a). However, it is still unclear whether such codes in humans indeed reflect the WM content.
Noninvasive human brain imaging research has initially been limited to the study of elevated activity during WM delay periods. However, delayed activity cannot dissociate between the representation of WM content and other content-independent functions that assist in retaining and manipulating mental representations (Baddeley, 2012; Riggall and Postle, 2012; Sreenivasan et al., 2014a; D'Esposito and Postle, 2015). Multivoxel pattern analysis (MVPA) for fMRI has recently received substantial attention in WM research because this approach allows patterns of local activity to be related directly to particular WM contents (Norman et al., 2006; Haynes, 2015; Lee and Baker, 2016; Christophel et al., 2017) and therefore to identify brain regions that code specific types of mental content. Controversially, these MVPA studies largely failed to find stimulus-specific activity in the PFC (Riley and Constantinidis, 2016). Instead, they suggest sensory and parietal regions to exhibit stimulus-specific activation (Lee and Baker, 2016), favoring sensory recruitment models of WM (Pasternak and Greenlee, 2005). Therefore, the view of the human PFC's function for WM has been reevaluated. It remains controversial if, or under which conditions, the PFC exhibits stimulus-specific codes (D'Esposito and Postle, 2015; Ester et al., 2015; Postle, 2015; Lee and Baker, 2016; Riley and Constantinidis, 2016).
Meanwhile, evidence has accumulated that different types of information are coded in different brain regions depending on their degree of abstraction (Christophel et al., 2017). Vibratory frequencies are not thought to be mere low-level sensory stimulus features but are considered as “analog properties” together with stimulus features such as intensity or duration. Most recent studies focused on sensory aspects of WM and mainly fostered sensory recruitment models of WM (Sreenivasan et al., 2014). Research on such more abstract (analog) features has been limited.
Based on the monkey work of Romo and colleagues (Romo and de Lafuente, 2013), we sought to investigate whether similar regions in the human brain also exhibit activation patterns to which MVPA is sensitive. To test which brain regions parametrically represent vibratory frequency, we applied a searchlight approach using support vector regression (Kriegeskorte et al., 2006; Kahnt et al., 2011) in a human version of the vibrotactile WM task. Crucially, this approach does not make any a priori assumptions on the localization of such representational codes.
Materials and Methods
Participants
All subjects (N = 24) were neurologically intact and right-handed, as assessed by the Edinburgh Handedness Inventory (EHI) (Oldfield, 1971). A total of N = 22 subjects were included in the analysis (age: 26.46 ± 3.38 years, 16 female, EHI: 0.87 ± 0.15) because two subjects were excluded due to excessive head motion (>9 mm). All participants provided written informed consent for the procedure in accordance with protocols approved by the local ethics committee of Freie Universität Berlin.
Experimental design
We applied a retro-cue WM paradigm in which the presentation of two vibrotactile sample frequencies was followed by a visual cue indicating which of the two frequencies had to be memorized during a subsequent delay period (Fig. 1). Subjects performed a two-alternative forced choice task after a delay of 12 s to indicate which of two test stimuli was identical to the memorized sample. Responses were given via right middle and index finger button presses and left/right response mapping was randomized across subjects. Vibrotactile stimuli were delivered to the left index finger using a 16-dot piezoelectric Braille-like display (4 × 4 matrix with 2.5 mm spacing) controlled by a programmable stimulator (QuaeroSys, St. Johann, Germany).
The set of vibrotactile sample stimuli contained four frequencies: 10, 22, 34, and 46 Hz (Fig. 1). All 16 pins of the Braille-like display vibrated with equal amplitude for 600 ms with smoothed onsets and offsets. One of the test stimuli matched the sample stimulus' frequency, whereas the foil stimulus was higher or lower in 50% of the trials, respectively. To achieve equal difficulty in each trial, we adjusted the frequency of the foil following the Weber–Fechner law (logarithmic relation between stimulus and perception; Fechner, 1966) anchored at 28 + 9 Hz, resulting in 7.6 and 13.2 Hz for the 10 Hz stimulus, 16.6 and 29.1 Hz for the 22 Hz stimulus, 25.7 and 44.9 Hz for the 34 Hz stimulus, and 34.8 and 60.8 Hz for the 46 Hz stimulus. The mask stimulus was composed of pins vibrating at all four sample frequencies and 180° phase-shifted equivalents and was presented for 500 ms with smoothed onsets and offsets.
Each experimental run consisted of 48 trials with a 12 s WM delay, supplemented by 8 catch trials with delays of 4 or 8 s to ensure sustained WM throughout the retention period. Each of six stimulus pairs was presented equally often in a balanced stimulus order, in which each of the four samples was memorized 12 times. Subjects were familiarized with the task by completing 1 or 2 practice runs outside of the scanner 1–5 d before the fMRI experiment.
fMRI data acquisition and preprocessing
fMRI data were acquired in 4 runs of 20 min on a 3 T scanner (TIM, TRIO, Siemens) at the MRI facility of the Freie Universtät, Berlin. For each run, 600 functional images were acquired (T2*-weighted gradient-echo EPI: 37 slices; interleaved order; 20% gap; whole brain; TR = 2000 ms; TE = 30 ms; 3 × 3 × 3 mm3 voxel; flip angle = 90°; 64 × 64 matrix). In addition, structural MRI data were acquired (MPRAGE, 176 sagittal slices, TR = 1900 ms, TE = 2.52 ms, 1 × 1 × 1 mm3 voxel).
Trial onsets were time locked to the functional images. This approach allows a TR-wise analysis of the 12 s WM phase in six consecutive functional images, as applied previously in the study of visual short-term memory (Christophel et al., 2012).
Statistical analysis
Time-resolved searchlight decoding of parametric multivariate frequency representations.
To preserve the spatiotemporal structure of the fMRI signal, no smoothing, normalization, or slice-time correction was performed before classification and preprocessing was limited to spatial realignment. fMRI data processing was performed using SPM8 (Wellcome Trust Centre for Neuroimaging, Institute for Neurology, University College London). Finite impulse response (FIR) models were used to obtain runwise beta estimates for each 2 s time bin of the WM delay (to display time courses, an additional 2 s time bin before and after the WM delay was modeled). High-pass filtered data (cutoff 192 s) were included in a model with a total of 132 beta estimates (4 stimuli × 8 time bins × 4 runs + 4 constants). All trials of the experiment except the catch trials were modeled to avoid any potential biases from the exclusion of trials. As control analysis, the identical decoding approach was performed with beta estimates from a model containing only correct trials and revealed virtually identical results.
To identify where in the brain information about the memorized vibratory frequency was encoded, we applied a time-resolved multivariate decoding analysis using a support vector regression (SVR) model, which has been described previously (Kahnt et al., 2011). This approach allows testing for local multivariate representations (activation patterns of voxels) that code the vibratory frequency maintained during WM. Compared with classic support vector machine approaches, SVR does not only distinguish between two classes of stimuli. Instead, an SVR model is trained to predict the value of a continuous variable, here vibratory frequency, based on a multivariate data vector. This is similar to a univariate regression approach, in which the prediction of a dependent variable is only based on one independent variable. In SVR, a data vector with multiple independent variables (activation pattern of voxels) is used to predict a single dependent variable (the maintained vibratory frequency).
All decoding analyses were performed using The Decoding Toolbox (TDT) (Hebart et al., 2015), which provides an interface to apply the computational routines of LIBSVM (Chang and Lin, 2011) and LIBLINEAR (Fan et al., 2008) to neuroimaging data. Crucially, we did not limit our analysis by assumptions on where in the brain local activation patterns relate to WM content. Instead, we applied a searchlight analysis (Kriegeskorte et al., 2006) in which a spherical volume of interest is moved step-by step (voxel-by-voxel) through the entire brain. In each searchlight step, a SVR model is trained and tested on data extracted from a local sphere of voxels. Thereby, 3-dimensional prediction accuracy maps for the whole brain are derived for each experimental subject. Consecutively, these maps allow group-level statistical testing to identify those regions that exhibit above-chance predictions across subjects and thus, maintain information on the memorized vibratory frequency.
The subject-level decoding analysis was performed using a leave-one-out cross-validation scheme in combination with an SVR kernel (linear kernel and constant regularization parameter of c = 1). Independent whole-brain searchlight decoding analyses were performed for each time bin. Beta estimates of each of the eight time bins were first z-scaled (normalized) across the samples for each voxel, as implemented in TDT. In each searchlight step, z-scaled beta estimates of voxels in a four-voxel radius sphere were extracted from each of the four maps corresponding to the four vibratory frequencies. Next, z-scaled beta estimates were forwarded to a fourfold cross-validation scheme in which data from three of the four runs are used to “train” the SVR. This SVR is then used to predict values of vibratory frequency based on the independent fourth run's data. To make our data comparable to previous reports, we used as a “test” measure the prediction accuracy, defined as the Fisher's z-transformed correlation coefficient between the predicted value levels and the actual value levels of the test dataset (Kahnt et al., 2011) (in TDT: “zcorr”). The center of the searchlight was moved voxelwise through the brain and prediction accuracy values were saved in the corresponding whole-brain accuracy maps. In this way, we gained a prediction accuracy map per subject reflecting which brain regions code the value of the vibratory frequency in a multivariate way, which is considered the multivariate pendant of a parametric code (Kahnt et al., 2011). This analysis yielded eight accuracy maps per subject, corresponding to the six timebins spanning the WM interval and one volume before and after it (Fig. 1).
Accuracy maps were normalized to MNI space using unified segmentation (as implemented in SPM8) and smoothed with a 3 mm full-width half-maximum kernel. These were entered consecutively into a flexible factorial second-level design (repeated measures) and a t-contrast for the mean of the six WM time bins against zero was computed. Voxels with above-chance prediction accuracies are reported at a threshold of p < 0.05 family-wise error (FWE) corrected for multiple comparisons and displayed with a cluster extent threshold of 35 voxels. All reported coordinates are in MNI space. The SPM anatomy toolbox was used to establish cytoarchitectonical references (Eickhoff et al., 2007).
Control analysis: decoding the nonmemorized stimulus.
As a control analysis, we tested for above-chance prediction accuracy for the nonmemorized stimulus. New FIR models were estimated, in which the four sets of FIR regressors modeled those trials in which a particular stimulus was presented but not memorized. Therefore, each beta image was estimated with equal amounts of data (modeling the same amount of trials) as in the original analysis. Beta images were entered into an identical SVR searchlight and second-level analysis as in the main analysis. No clusters with above-chance prediction accuracies were found.
Control analysis: label permutation tests.
To determine the prediction accuracy of an SVR analysis, the beta estimates within a searchlight volume for all four ordered stimulus conditions are entered into a SVR model. A high prediction accuracy is only expected if the activation patterns in a given region represent this order of the four frequencies. To test for the specificity of our analysis, we performed exhaustive whole-brain SVR searchlight analyses for all possible permutations of frequency labels. To this end, we determined for all possible permutations the distance from the rank order calculated as the sum of absolute difference of adjacent ranks (e.g., the linear order of 10, 22, 34, 46 Hz has a distance of ranks of sum(|1 − 2| + |2 − 3| + |3 − 4|) = 3 and the permutation 22, 10, 34, 46 Hz corresponds to the sum(|2 − 1| + |1 − 3| + |3 − 4|) = 4, resulting in a difference of 1 from the linear order). Because rank order is symmetric, the permutation analyses span 12 instead of 24 possible permutations congregated into four classes of distances from the linear order. We extracted the prediction accuracies from the group peak voxels of the original analysis. Figure 2C presents averaged prediction accuracies for the four groups of distances from the linear order.
For a statistical assessment, we tested for each time bin whether an increase in ordering correlated with an increase in prediction accuracy. We calculated Spearman correlation coefficients for each time bin and report significant effects at p < 0.05 (Figure 2C).
Results
Behavioral performance
Subjects (N = 22) successfully memorized the relevant information in the demanding task, as indicated by the 62 ± 4% (mean ± SD) correct responses over 4 fMRI runs. As expected, performance was better for the second (64.5 ± 5%) than for the first (60 ± 7%) presented stimulus (paired t-test t(21) = −2.65, p = 0.015), which was accounted for by a counterbalanced stimulus presentation.
Multivariate mapping of regions that code the content of WM
To test which brain regions exhibit a WM code for vibratory frequency, we conducted a time-resolved whole-brain SVR analysis. We applied an SVR approach to test in which regions information on the order of increasing frequencies can be found by testing where in the brain decoding performance exceeded chance level. We computed a t-contrast of prediction accuracy versus zero across the WM delay period within a second-level flexible factorial design. After accounting for multiple comparisons with a conservative thresholding of p < 0.05 and FWE correction, we identified five clusters of above chance decoding that match with reports from monkey research: bilateral dorsal premotor cortex (dPMC), a cluster in the supplementary motor area/cingulate cortex (SMA/CC), and a cluster in the right inferior frontal gyrus (IFG), as displayed in Figure 2A and Table 1.
As a control analysis, the study design allowed us to perform an identical MVPA for the noncued stimulus because two sample stimuli were presented on each trial and pairs of stimuli were balanced (Christophel et al., 2012). Decoding of the nonmemorized stimulus did not reveal any significant clusters of voxels, demonstrating the specificity of the main analysis.
To test for the specificity of the linear ordering of stimuli with the applied SVR approach, we performed permutation tests. To this end, we relabeled the training and test data with all possible permutations of frequency labels. If the order of frequencies were crucial for the SVR's performance, then one would expect increased prediction accuracies with increasing order. We particularly expected this effect during the WM delay period, but not before or after it. The time courses for all four types of permutations across the delay period are plotted in Figure 2C. Note that completely unordered labeling leads to chance performance of the SVR throughout the entire WM period. With increased order, the prediction accuracy increased for all identified regions across the delay period. As expected from the dynamics of the BOLD response, we found that this effect peaked between 4 and 8 s after the WM onset.
In addition to the multivariate analysis, we performed a classic general linear model analysis and tested for univariate parametric modulations during the WM delay phase, which did not reveal any significant voxels (p < 0.001, uncorrected). This suggests that it is the combination of multiple voxels, and thus, distributed neuronal populations, that code for the vibratory frequency rather than an overall parametrically increased activity throughout voxels in the identified clusters.
Discussion
Studies in nonhuman primates (Romo et al., 2012) have identified a network of brain regions where neurons exhibit firing behaviors that are parametrically modulated by the maintained vibrotactile frequency. In a series of EEG experiments, Spitzer and coworkers extended this work to humans by showing that, consistent with the monkey reports, the retention of vibration frequency induces a content-specific parametric beta-band modulation during delay periods, which was source-localized to the right PFC (Spitzer and Blankenburg, 2012; Spitzer et al., 2014a). However, until now, the mapping of brain regions that code vibrotactile frequency in the form of multivariate patterns during human WM has been lacking. Here, we used a multivariate regression approach to test which brain regions exhibit specific patterns of activation for maintained vibratory frequencies. The applied SVR approach can be considered as the multivariate equivalent of univariate tests for parametric modulations; that is, a classifier is trained to predict a value (here: the frequency) instead of a class label (Fan et al., 2008; Kahnt et al., 2011; Hebart et al., 2015). Therefore, the order of frequencies needs to be represented in a region to yield above-chance decoding results. With a whole-brain searchlight approach, testing within the entire brain for such local content-specific codes, we found that the bilateral premotor cortices, medial frontal regions, and right IFG code frequency-specific information. Our control analyses proved the specificity of these results because it was not possible to detect information on the nonmemorized stimuli and because permutation tests showed that it was indeed the ordering of the frequencies that enabled to decode information from these regions.
Interestingly, the identified network overlaps with brain regions reported in Romo and colleagues' monkey research (Romo and de Lafuente, 2013; Romo et al., 2012), as well as with those in the human EEG reports (Spitzer and Blankenburg, 2012; Spitzer et al., 2014a). Here, we extend these findings by showing for the first time that the right PFC also codes multivariate content-specific stimulus information.
In parallel to the monkey data, WM content was not limited to a single brain region; rather, a network of regions was identified. It is reasonable to assume that the brain represents mental content depending on task demands in distributed cortical networks, which can also exhibit different types of codes (Schlegel et al., 2013; Larocque et al., 2014; Lewis-Peacock et al., 2015; Postle, 2015; Lee and Baker, 2016; Christophel et al., 2017). In addition to the PFC, our analysis revealed that the bilateral dorsal premotor cortices and supplementary motor regions (Romo and Salinas, 2003; Hernández et al., 2010; Romo and de Lafuente, 2013; Romo et al., 2012) also exhibit these codes. Dorsal premotor cortices have been indicated to contribute to rehearsal processes (Fegen et al., 2015), whereas details on the mechanisms underlying rehearsal or refreshment of sensory contents (Chun and Johnson, 2011) are still relatively unclear (Baddeley, 2012). Our study shows that these regions exhibit a parametric representation of the WM content. Even though this is consistent with the monkey reports, a direct mapping between monkey and human anatomy remains speculative until future studies provide more direct evidence for these regions' functional roles. Although our study provides a first link for the apparent gap between unit-level recordings and fMRI-based classification approaches, MVPA studies are limited in their ability to identify a specific coding schema. They also allow us to test which regions code WM information but cannot resolve the nature of the code as such.
Our findings are in contrast to several recent WM MVPA studies. Most of these studies have argued for sensory recruitment models because they found content-specific codes predominantly in sensory and posterior regions (Pasternak and Greenlee, 2005; Sreenivasan et al., 2014; Riley and Constantinidis, 2016). Similarly, early monkey research on tactile WM supported the notion that sensory regions are predominantly recruited to represent WM content (Zhou and Fuster, 1996, 2000). This view aligns with the famous WM model by Alan Baddeley (Baddeley, 2012), which is based on the idea of a visuospatial sketchpad, which is now often equated with sensory regions that are supposed to act as a kind of whiteboard on which mental content is drawn (Pasternak and Greenlee, 2005). In addition to WM studies, research on mental imagery also builds heavily on the idea that activity in sensory regions represents the content of mental images (Kosslyn, 2005; Tong, 2013). Indeed, it has been demonstrated that primary somatosensory cortex is recruited for tactile imagery (Schmidt et al., 2014).
On the contrary, frontal and prefrontal regions are well documented to activate in WM studies across different modalities, including tactile WM (Preuschhof et al., 2006). However, these regions are most often considered to fulfill supporting, top-down functions rather than to represent WM content (D'Esposito and Postle, 2015; Myers et al., 2017). With the exception of studies on vibrotactile frequency representations, only a few fMRI studies have revealed prefrontal WM codes that reflect aspects of the WM content (Lee et al., 2013; Ester et al., 2015).
One interpretation of these conflicting findings might be that different tasks require various types of codes. Visual WM studies mostly necessitate mentally representing very precise visuospatial stimulus features such as orientation (Albers et al., 2013), position (Jerde et al., 2012), motion direction (Emrich et al., 2013; Christophel and Haynes, 2014), or the layout of a stimulus (Christophel et al., 2012; Lee et al., 2013). These studies convincingly demonstrate that such low-level or sensory-like features are maintained in brain regions which also deal with such type of information during perception. Therefore, Christophel and coworkers recently suggested that it might be the level of “abstraction” of the WM content that determines the network of regions coding it (Christophel et al., 2017). Although primary sensory cortices deal with low-level sensory-like material, regions higher up in the cortical hierarchy, such as the PFC, might process more abstract stimulus features. This view resonates well with evidence that the right PFC processes nonverbal quantitative types of information (Nieder and Dehaene, 2009) and the left lateral PFC, particularly BA44, is well known for its pivotal role in language processing. In monkeys, the right PFC has also been associated with codes for size/length/quantity or numerosity (Nieder and Dehaene, 2009; Nieder, 2013). These magnitude-related stimulus properties are considered to be represented independently of sensory aspects (Walsh, 2003) and human fMRI and EEG studies also associate them with the PFC (Pinel et al., 2004; Nieder and Dehaene, 2009; Spitzer et al., 2014b, 2014c). Consistent with this reasoning, the maintenance of vibrotactile frequency might be considered a representation of an abstract, magnitude-like representation despite a low-level stimulus feature, a type of code for hierarchically higher brain regions such as the PFC. Even though this idea is appealing and in accordance with our data, future research has to provide evidence for it by dissecting different types of information and codes, as well as their functional contribution to goal-directed purposeful behavior and decision making (Romo et al., 2012), including other modalities.
Conclusion
Our study found the right PFC, SMA/CC, and bilateral dPMC to contain multivariate codes of retained vibratory frequencies during human WM. These codes represented the parametric order of the memorized frequencies and the identified network overlapped with previous monkey studies. Future studies will elicit what types of WM contents are represented differentially across species and which cellular mechanisms can be translated more directly between them. For this endeavor, it is crucial to establish links between different observational levels in neuroscience. The bridges between monkey and human data in the study of WM are rare. Our results allow initial new perspectives to unite the different lines of literature using distinct methodologies.
Footnotes
T.T.S. was supported by the Research Training Group (Grant GRK 1589/2, Sensory Computation in Neural Systems, by the German Research Foundation, DFG).
The authors declare no competing financial interests.
- Correspondence should be addressed to Timo Torsten Schmidt, Department of Education and Psychology, Neurocomputation and Neuroimaging Unit, Freie Universität Berlin, Habelschwerdter Allee 45, 14195 Berlin, Germany. timoschm{at}uos.de