Parametric Representation of Tactile Numerosity in Working Memory

Abstract Estimated numerosity perception is processed in an approximate number system (ANS) that resembles the perception of a continuous magnitude. The ANS consists of a right lateralized frontoparietal network comprising the lateral prefrontal cortex (LPFC) and the intraparietal sulcus. Although the ANS has been extensively investigated, only a few studies have focused on the mental representation of retained numerosity estimates. Specifically, the underlying mechanisms of estimated numerosity working memory (WM) is unclear. Besides numerosities, as another form of abstract quantity, vibrotactile WM studies provide initial evidence that the right LPFC takes a central role in maintaining magnitudes. In the present fMRI multivariate pattern analysis study, we designed a delayed match-to-numerosity paradigm to test what brain regions retain approximate numerosity memoranda. In line with parametric WM results, our study found numerosity-specific WM representations in the right LPFC as well as in the supplementary motor area and the left premotor cortex extending into the superior frontal gyrus, thus bridging the gap in abstract quantity WM literature.


Introduction
Humans can tell whether 100 people are a larger group than 50 people quite accurately without counting. This ability to quantify amount, size, length, or other analog stimulus properties can be performed nonsymbolically, independent of language (Dehaene, 1992;Spitzer et al., 2014b). Indeed, human infants and several animals are able to approximate a variety of quantities (Nieder, 2005;Piazza et al., 2007;Piazza and Izard, 2009;Nieder and Dehaene, 2009), suggesting a common elemental system which has been termed the approximate number system (ANS; Gallistel and Gelman, 1992;Dehaene, 2011).
While numerosity is a discrete stimulus property, the ANS allows an approximation of numerosity, resulting in an analog estimation. Thus, in contrast to the symbolic mental representation of numbers as classes or categories, it has been hypothesized that the ANS representation resembles that of continuous quantities or magnitudes such as intensities, lengths, or frequencies (Piazza et al., 2004;Nieder and Dehaene, 2009;Spitzer et al., 2014a). In support of this, neural representations underlying both the ANS and continuous quantities have been shown to be supramodal, implying a representation abstract in nature (Piazza et al., 2006;Spitzer and Blankenburg, 2012;Spitzer et al., 2014a;Vergara et al., 2016). Moreover, small numbers are rapidly and accurately identified without counting, known as subitizing (Kaufman et al., 1949). Thus, these numbers are represented as discrete values. If the number of items exceeds the subitizing threshold, counting is required to determine the exact amount. When there is insufficient time for counting, the ANS approximates the quantity in a fast and efficient manner.
The functional anatomy of the ANS has been extensively characterized in both human and nonhuman primates (NHPs). A frontoparietal network comprising the dorsolateral prefrontal cortex and the posterior parietal cortex (PPC), specifically the intraparietal sulcus (IPS), is involved in approximating quantities during perception Piazza et al., 2004Piazza et al., , 2007Cantlon et al., 2006Cantlon et al., , 2009Jacob and Nieder, 2009;Knops and Willmes, 2014). Moreover, the right hemisphere has been shown to be dominant with respect to quantity estimation (McGlone and Davidson, 1973;Young and Bion, 1979;Kosslyn et al., 1989); however, recent studies have found that both hemispheres respond to approximate visual numerosity (Piazza et al., 2004;Ansari et al., 2006). Particularly in nonsymbolic numerosity perception, the IPS has been shown to exhibit stronger numerosity-selective responses than the PFC (Tudusciuc and Nieder, 2009), and the PPC, especially the IPS, responds to the nonsymbolic numerosity processing (Piazza et al., 2004(Piazza et al., , 2007. The ANS literature is primarily focused on perception with relatively few NHP studies extending to investigate working memory (WM) representations of approximate quantities (Nieder, 2016). As short-term maintenance of information is critical for higher-order cognitive functions such as decision-making and reasoning, it is crucial to investigate beyond perception to the maintenance of approximate quantities in WM. In line with results from perception studies of the ANS, neurons in the frontoparietal network were found, specifically in the PFC and IPS, to exhibit numerosity-selective activity during WM (Jacob et al., 2018). Furthermore, supramodal coding of numer-osity memoranda in the frontoparietal cortex has been identified (Nieder, 2017). Interestingly, in contrast to perception, the proportion of numerosity-selective neurons in the PFC and their tuning strength to numerosity have been more prominent than the ones in the PPC during WM retention. Moreover, neurons in the PFC remained selective and discriminated numerosities better than neurons in the PPC during the WM delay (Nieder and Miller, 2004;Tudusciuc and Nieder, 2009;Nieder, 2016).
To the best of our knowledge, only a single study has focused on the WM representation of numerosity in humans (Spitzer et al., 2014a), although some approximate numerosity perception studies used fMRI multivariate pattern analysis (MVPA) method with WM-related paradigms focusing on the perceptual processes instead of the WM retention (Eger et al., 2009;Borghesani et al., 2019;Castaldi et al., 2019). Spitzer et al. (2014a) probed the oscillations underlying multimodal WM representations by training participants to estimate numerosity from sequential auditory, visual, and tactile stimuli. They identified strong and long-lasting alpha oscillations in the PPC reflecting WM load, whereas, in line with NHP results, beta-band activity in the right PFC showed numerosityselective modulation.
Nevertheless, whole-brain research regarding the localization of numerosity memoranda in humans is lacking. To this end, we designed a tactile delayed match-tonumerosity (DMTN) task in combination with whole-brain, searchlight, MVPA of human fMRI data (Christophel et al., 2012;Schmidt et al., 2017;Uluç et al., 2018). Using this analysis approach, we localized brain regions maintaining approximate number content in WM. As per previous studies (Spitzer et al., 2014a;Nieder, 2016), we hypothesized that the content would be represented in frontal regions, specifically the right PFC.

Participants
Thirty-eight healthy volunteers participated in the study. The sample size was based on the successful use of similar sample sizes in earlier MVPA experiments with analog experimental designs and analyses (Schmidt et al., 2017;Christophel et al., 2018). In addition, it accords with recent theoretical work on power analysis for random field theory-based cluster-level statistical inference (Ostwald et al., 2019). The data of four participants were excluded due to low performance levels (Յ60%), resulting in data from 34 participants (mean Ϯ SD age, 25.53 Ϯ 5.43 years; 19 females) being further analyzed. All were right handed according to the Edinburgh Handedness Inventory with a mean Ϯ SD index of 0.82 Ϯ 0.14 (Oldfield, 1971). The experimental procedure was approved by the local ethics committee and was conducted in accordance with the Human Subject Guidelines of the Declaration of Helsinki. All participants provided written informed consent before the experiment and were compensated for their participation.

Stimuli
Tactile stimuli consisted of trains of square-wave electric pulses (200 s) delivered via a pair of surfaceadhesive electrodes attached to the participant's left wrist. A constant current neurostimulator (model DS7A, Digitimer) was used to deliver the stimuli. Subjects reported tactile sensations radiating to the thumb, index, and middle finger, verifying stimulation of the median nerve. Individual sensory thresholds were determined for each participant. The stimulus intensity was then adjusted to a target value of ϳ200% of the sensory threshold (mean, 6.42 mA; SD, 1.20 mA).
A to-be-remembered stimulus sequence comprised 7, 9, 11, or 13 pulses. To dissociate stimulus length and perceived pulse frequency (spacing of tactile pulses) from the numerosity of pulses, the duration of the stimulus varied, and the interpulse intervals were randomized. To this end, we defined four stimulus durations (960, 1020, 1080, and 1140 ms). Each duration was subdivided into 60 ms slots, resulting in 17, 18, 19, and 20 slots, respectively. The temporal distribution of the pulses was then randomized across the slots (Fig. 1A, illustrative stimuli). Within each run, each numerosity was presented in a short (17 or 18) and a long (19 or 20) duration, resulting in 24 different numerosity-duration pairings (4 numerosities ϫ 2 durations/run ϫ 3 uncued numerosities). The different durations were balanced across runs. The alternatives for each cued numerosity were computed according to the respective sample (Ϯ3 pulses). Additionally, the target stimulus and the cued sample never had the same duration, ensuring that memorizing the duration or average frequency of the target does not help to perform the task. We also performed a Fourier transformation of the stimuli, which ensured that all stimuli were composed of similar combinations of frequencies. Therefore, this stimulus design ensured that participants had to memorize the stimulus numerosity since they could not use the temporal density of the pulses or the stimulus length as WM memoranda to solve the task.

Task
We used a DMTN paradigm in which participants remembered the estimated numerosity of a stimulus. Each trial began with the presentation of two pulse sequences with different numerosities. Next, a retro-cue ("1" or "2") indicated which of the two numerosities had to be remem- Figure 1. Sample pulse sequences and experimental paradigm. A, Sample stimuli. Pulse sequences of 7, 9, 11, and 13 were used as experimental stimuli. For each numerosity, there were four different durations (960, 1020, 1080, and 1140 ms), where each duration was subdivided into 60 ms slots. The distribution of pulses to slots was randomized for each stimulus presentation. The first and the last slot of each stimulus always contained a pulse. The stimuli displayed are for illustrative purposes. B, Experimental paradigm. A delayed match-to-numerosity task was used, where two sample stimuli and a mask were presented consecutively. A visual retro-cue that was presented simultaneously with the mask indicated which of the numerosities should be retained for the 12 s delay. After the delay, participants performed a two-alternative forced choice, indicating which of the two test stimuli had the same numerosity as the cued stimulus. The response period was 1.5 s. Please note that the stimulus duration and interstimulus interval changed depending on the stimulus duration, but the onset of each event was locked to coincide with the onset of an image acquisition. bered. To suppress potential perceptual residues, in the sense of afterimages (Sperling, 1960;Christophel and Haynes, 2014;Christophel et al., 2015), a mask consisting of the longest duration (1140 ms) with a pulse in each of the 20 slots, was applied simultaneously with the onset of the retro-cue. Following a 12 s retention phase, two test stimuli were presented and a two-alternative forced choice was given. Neither of the test stimuli were identical to the encoded stimulus; however, one had the same numerosity, while the duration and the frequency were different. This ensured that participants used the approximated numerosity of the stimulus instead of some other stimulus feature to correctly match the test with the remembered stimulus. The numerosity of the alternative stimulus was three pulses plus or minus the target stimulus. To ensure that the number of pulses in a sequence could not be easily counted, the lower alternative stimulus for the lowest to-be-remembered numerosity (7), was set to 5 and thus above a previously established subitizing threshold of approximately 4 (for tactile modality, it was shown to be 3-4; Riggs et al., 2006;Plaisier et al., 2009Plaisier et al., , 2010Plaisier and Smeets, 2011;Spitzer et al., 2014a;Tian and Chen, 2018). After the second target stimulus, participants had 1.5 s to indicate, via button press with their right middle or index finger, which of the two stimuli had the same numerosity as the encoded stimulus ( Fig.  1B, experimental design). Furthermore, the response mapping was counterbalanced across participants. In total, a trial lasted 21 s and an experimental run, consisting of all possible stimulus pairings presented equally often (12 pairings ϫ 4 presentations ϭ 48 trials) in a randomized order, with intertrial intervals of 1.5 or 3.5 s, lasted 18.7 min. Four experimental runs were collected for each participant, resulting in a total recording time of 74.8 min.
Before the fMRI experiment, each participant was familiarized with the timing and structure of the task by performing up to two experimental runs outside the scanner.

Number naming test assessing countability
Subsequent to the fMRI session, we applied a numbernaming task to ensure that participants were unable to count the number of pulses used in the stimulus set.
Participants were asked to try to count the number of pulses. The stimuli ranged from 1 to 15 pulses with 5 different duration and temporal pulse distribution combinations of each numerosity were tested, resulting in 75 trials. The counting test was performed after fMRI data acquisition so as to prevent biasing the participants toward counting the pulses in the main experiment.
To ensure that the presented numerosities were above participants' subitizing thresholds, we calculated the mean performance for each numerosity across participants and calculated each average estimated numerosity. We then compared the slope of accuracy for estimating numerosities with earlier studies that calculated subitizing thresholds for tactile stimuli (Riggs et al., 2006;Plaisier et al., 2009Plaisier et al., , 2010Plaisier and Smeets, 2011;Spitzer et al., 2014a;Tian and Chen, 2018). We performed a linear trend analysis using linear regression to determine whether the distance between the true and estimated numerosity scales with increasing true numerosity in a linear fashion.

fMRI data acquisition and preprocessing
fMRI data were acquired in four runs, with a Siemens 3 T Tim Trio MRI scanner (Siemens) equipped with a 32channel head coil. In each run, 565 images were collected (T2‫-ء‬weighted gradient echo EPI: 37 slices; ascending order; 20% gap; whole brain; TR ϭ 2000 ms; TE ϭ 30 ms; 3 ϫ 3 ϫ 3 mm³; flip angle ϭ 70°; 64 ϫ 64 matrix). After the last functional run, a high-resolution structural scan was recorded using a T1-weighted MPRAGE sequence (1 ϫ 1 ϫ 1 mm³; TR ϭ 1900 ms; TE ϭ 2.52 ms; 176 sagittal slices). fMRI data preprocessing was performed using SPM12 (Wellcome Trust Center for Neuroimaging, Institute for Neurology, University College London, London, UK). Functional images were slice time corrected and spatially realigned to the mean image. To conserve the spatiotemporal structure of the fMRI data for the multivariate analyses, no smoothing or normalization was performed. For the univariate control analysis, functional images were normalized to MNI space and smoothed with an 8 mm FWHM kernel.

First-level finite impulse response models
A time-resolved, multivariate searchlight analysis (Kriegeskorte et al., 2006;Schmidt et al., 2017) was used to identify brain regions encoding memorized numerosity information. First, a general linear model (GLM) with a set of finite impulse response (FIR) regressors was fit to each participant's data to obtain runwise parameter estimates of each WM content (numerosity value of 7, 9, 11, or 13). A single FIR regressor was estimated for each fMRI image or 2 s time bin (1 TR); thus, the 20 s trial was divided into 10 time bins. We additionally included the first five principal components accounting for the most variance in the CSF and white matter signal time courses, respectively (Behzadi et al., 2007), and six head motion regressors, as regressors of no interest. Moreover, the data were filtered with a high-pass filter of 128 s. The resulting parameter estimates were used for the MVPA, performed with The Decoding Toolbox (TDT) version 3.52 .

Multivariate pattern analysis
For the decoding of memorized numerosity information, a searchlight-based multivariate analysis using a support vector regression (SVR) approach was performed with the computational routines of LIBSVM (Chang and Lin, 2011), as implemented in TDT. SVR MVPA (for more discussion, see Kahnt et al., 2011;Schmidt et al., 2017) considers the variable of interest (memorized numerosity) as a continuous data vector with multiple independent variables (multivariate BOLD activities) as opposed to the commonly used support vector machine approach that treats the variable of interest as a categorical object. This means that the SVR MVPA approach seeks a linear continuum for the numerosities in which their distance is proportional to the distances of the rank order.
We analyzed each time bin independently by implementing a searchlight decoding analysis with a spherical searchlight radius of 4 voxels. For a given voxel, z-scaled parameter estimates (across samples) corresponding to each WM condition were extracted from all voxels within the spherical searchlight for each run. This yielded 16 pattern vectors (4 WM contents ϫ 4 runs), each corresponding to the BOLD activity pattern for a specific WM condition of a functional run. We then fitted a linear function to these pattern vectors such that the multivariate distribution for the different numerosities follows a linear mapping of numerosities. The z-scaled parameter estimates were entered into an SVR model with a fixed regularization parameter c that was set to 1.
We used a leave-one-run-out cross-validation scheme for the subject-level decoding analysis. The SVR classifier was trained on three runs (12 pattern vectors) and tested on the data of the independent fourth run (4 pattern vectors) for how well it predicted the values of the remaining run. The allocation of training and test runs was iterated so that each of the four functional runs was used as a test run once, resulting in four cross-validation folds. The prediction performance from each cross-validation fold was reported by a Fisher's z-transformed correlation coefficient between the predicted and the actual numerosity information estimate. The mean prediction accuracy across cross-validation folds was assigned to the center voxel of the searchlight, and the center of the searchlight was moved voxel by voxel through the brain, resulting in a whole-brain prediction accuracy map. Consequently, we obtained one prediction accuracy map for each time bin for each participant, where the prediction accuracy reflects how well a linear ordering according to the associated numerosities could be read out from the locally distributed BOLD activity pattern at a given voxel location and time.
Next, prediction accuracy maps were normalized to MNI space and smoothed with an 8 mm FWHM kernel. They were then entered into a second-level, repeatedmeasures ANOVA with subject and time (time bins) as factors. To assess which brain regions exhibit WM content-specific activation patterns during the delay period, we computed a t-contrast across the six time bins corresponding to the 12 s WM delay (time bins 3-8). The results are presented at p Ͻ 0.05 family-wise error (FWE) correction at the cluster level with a cluster-defining threshold of p Ͻ 0.001. Cytoarchitectonic references are based on the SPM anatomy toolbox where possible (Eickhoff et al., 2005). Presented images (e.g., surface projections with applied color scales) were created using MRIcron version9/9/2016 (McCausland Center for Brain Imaging, University of South Carolina, Columbia, SC).

Control analyses
In the first control analysis, we examined whether the decoded numerosity information during WM retention was specific to WM or could be assigned to perceptual residues. To this aim, we defined a second, first-level model with FIR regressors for the nonmemorized stimulus. We then implemented the identical searchlight decoding procedure as the main analysis. Thus, this control analysis tested for the presence of numerosity information of the nonmemorized stimulus.
Next, we conducted a parametric univariate analysis to ensure that the decoded information in the main analysis is not due to the modulation of mean activity level. To this end, we fitted a standard GLM with the following four HRF-convolved regressors: one regressor to capture WM processes, a parametrically modulated regressor for the numerosity content of the WM memoranda as well as eight [4 numerosities ϫ 2 (sample, test)] additional parametrically modulated regressors for each sample and test stimulus. First-level baseline contrasts for the parametric effect of memorized numerosity were forwarded to a second-level one-sample t test.
Finally, to test the specificity of the SVR analysis to the parametric order of the four numerosities, we performed exhaustive whole-brain SVR searchlight analyses for all possible permutations of numerosity labels. To achieve this, we computed distance rank order as a sum of the absolute difference of adjacent ranks [e.g., 11, 13, 7, and 9 numerosity is distance 5 (|3-4|ϩ|4 -1|ϩ|1-2|)] for all possible permutations of the numerosity order. Then, the permutations were grouped according to their distance from the original rank order. We used 12 instead of 24 permutations as the distances of rank order permutations are symmetric. Including the permutation with the correct linear order, the 12 permutations are aggregated into five classes depending on their distance from the correct linear order. Then, for each permutation analysis, we extracted the prediction accuracies of the group-peak voxels that are defined in the original analysis. For statistical assessment, we calculated the mean prediction accuracy across related time bins (WM time bins 3-8) for each peak voxel for each distance group (see Fig. 3C).

Behavioral performance
Thirty-four participants performed with 65.36 Ϯ 3.29% (mean Ϯ SD) accuracy in the demanding DMTN task across the four experimental runs ( Fig. 2A). To test whether the behavioral performance differed for the four numerosity values, we performed a one-way repeatedmeasures ANOVA with four levels, one for each numerosity. This test revealed a significant main effect (F (3,135) ϭ 7.52, p Ͻ 0.001). Post hoc t tests (Bonferroni corrected for multiple comparisons) between performances were significant for numerosity values 7 and 13 and 9 and 13 (p Ͻ 0.05/6; Fig. 2A). This is expected because we did not control for the Weber-Fechner effect except for the lowest numerosity (which we did due to subitizing concerns). As a result, as the numerosity increases, it becomes more difficult to differentiate between the sample and alternative stimuli, thus resulting in a lower performance for high numerosities (Fechner, 1966) but is unlikely to affect WM processing.

Behavioral performance on number naming test assessing countability
To test whether participants were able to count the numerosities used in the current study, participants per-formed an additional number-naming test. Previous research in tactile numerosity indicated the subitizing threshold for comparable stimuli to be four pulses (Riggs et al., 2006;Plaisier et al., 2009;Plaisier and Smeets, 2011;Spitzer et al., 2014a;Tian and Chen, 2018). The approximation of the subitizing threshold identified in the present study is in line with these reports (Fig. 2B). As expected, participants' perceptual accuracy decreased with increasing numerosity, and performance decreased to 50% when more than three pulses were presented. Similarly, the distance between the true and estimated numerosity increased with increasing numerosities (p Ͻ 0.001, linear trend analysis; Fig. 2C).

Multivariate mapping of regions that code numerosity as WM content
The time-resolved, searchlight-based multivariate regression analysis was performed to identify brain regions representing estimated numerosity memoranda. The SVR MVPA analysis for the WM retention period revealed numerosity-specific responses in the left premotor cortex (PMC) slightly extending to the primary motor area (MI), left middle frontal gyrus (MFG), left superior frontal gyrus (SFG) extending into bilateral supplementary motor areas (SMA), right SFG extending to the right frontal pole, and right MFG extending into the pars triangularis of the right inferior frontal gyrus (IFG). Results are reported at p Ͻ 0.05, FWE corrected at the cluster level with a clusterdefining threshold of p Ͻ 0.001 (Fig. 3, Table 1).
For the sake of completeness, we investigated whether numerosity information could be decoded from the IPS at an uncorrected statistical threshold of p Ͻ 0.001. We found a cluster in the right PPC extending to the IPS (peak at MNI: x ϭ 36, y ϭ Ϫ52, z ϭ 36 mm; z score ϭ 3.89; k ϭ 164), which was identified as hIP1 with a 39.5% probability and hIP3 with a 5.9% probability using the SPM anatomy toolbox (Eickhoff et al., 2005) at p uncorrected Ͻ 0.001.

Control analyses
To test whether the identified decoded information is indeed specific to the memorized numerosity representation, we applied the same searchlight procedure to the nonmemorized numerosity stimulus. This analysis did not reveal any clusters with above-chance prediction accuracy at p FWE-Cluster Ͻ 0.05.
Additionally, we conducted a univariate parametric analysis to test whether the decoding results could be due to differences in activation strength between WM contents. A second level t test revealed no significant voxels at p FWE-Cluster Ͻ 0.05, thus providing evidence for the multivariate nature of the numerosity representations identified in this study rather than the modulation of univariate mean activity.
Finally, we performed label-permutation tests to ensure the specificity of the linear ordering of stimuli in the SVR MVPA. Higher prediction accuracies were expected when the activation patterns in a given brain region represented the correct order of the four numerosity labels, and it was expected to decrease with the distance from the correct ordering. As expected, the prediction accuracy during WM was the highest for the true-labeled data and decreased with increasing distance from the correct ordering (Fig. 3C).

Discussion
The current study, to our knowledge, is the first to identify brain regions that code approximate numerosity WM content using human neuroimaging methods. Thus, this study extends the broad literature on ANS perception to the maintenance of mental representations, which can be used for higher-order cognitive functions. We used a well established, whole-brain, searchlight, DMTN paradigm to identify representations of tactile approximate numerosity memoranda. Specifically, we used an SVR technique, which, in contrast to support vector machines, treats the retained WM content as a continuous variable and thus predicts the ordering of content along the variable, rather than a singularly specific class label. Consequently, an above-chance prediction accuracy in a brain region means that the content-specific activation patterns follow a linear ordering according to the associated numerosity. Our searchlight analysis identified a distributed network spanning the left PMC, bilateral SFG, bilateral SMA, and right MFG extending into right IFG. Therefore, Figure 2. A, Mean rate of correct responses across participants (n ϭ 34) for different numerosities in the main WM DMTN task. The figure shows that the WM performance decreases with increasing numerosity. Error bars represent standard derivation (SD). Asterisks indicate statistical significance for pairwise t tests, Bonferroni corrected for multiple comparisons (p Ͻ 0.05/6). B, Mean performance across subjects for estimated numerosity in number naming task (mean Ϯ SD). C, True numerosities versus mean numerosity estimations (error bars show SD). these regions contain linearly ordered, multivariate WM representations of the numerosities.
Our results are in line with previous numerosity WM studies in NHPs and human EEG, which have established the central role of the PFC. Indeed, previous unimodal and multimodal studies have identified content-specific representations in the PFC (Nieder and Miller, 2004;Tudusciuc and Nieder, 2009;Spitzer et al., 2014a;Nieder, 2016;Jacob et al., 2018). More specifically, in humans, parametric modulation of upper-␤ oscillations in the right lat-eral PFC has been shown to reflect analog numerosity estimation that has been derived from discrete sequences, both within and between stimulus modalities (Spitzer et al., 2014a). Thus, the numerosity representations in the PFC are likely to be supramodal in nature. However, those studies used either electrophysiological recordings from an a priori brain region or have used univariate data analysis methods. The present study extends the literature on numerosity WM in the following two ways: first, to whole-brain fMRI data; and second, to Figure 3. A, Brain regions coding information for the memorized estimated numerosities. Group-level results of a t-contrast testing the 12 s WM delay for above-chance prediction accuracy. Brain regions carrying information about memorized scalar magnitudes are as follows: IFG, MFG, PMC, SMA, and SFG. B, Time courses of decoding accuracies of remembered (red) and nonremembered (gray) stimuli for all identified brain regions in the main analysis (Fig. 3A). Error bars indicate standard error of the mean (SEM). The figure shows that, for all clusters depicted in the main analysis, there is more numerosity-specific WM information for the remembered than for the forgotten numerosity, and the information is present throughout the WM delay period. C, Results of the label permutation tests. Five bars are shown for each brain region, respectively. Each bar displays the mean prediction accuracy estimated from the distance to correct order groups. The shade of the bar color, ranging from black to white, depicts the different distance to correct ordering. Black bars indicate the mean prediction performance of the group with the correct linear order, while white bars represent the mean prediction accuracy derived from the most linearly unordered data. Brain regions tested for label permutation are: IFG, MFG, PMC, SMA, and SFG. Error bars indicate SEM. multivariate data analysis methods, specifically the SVR MVPA. The benefits of multivariate over univariate analysis methods have been well established (Haynes, 2015). Multivariate analysis techniques are sensitive to the combinatorial aspects of voxel activity, thereby enabling the identification of spatially distributed representations (Haynes, 2015;Hebart and Baker, 2018). Thus, our results agree with and extend the previous NHP and human EEG numerosity WM findings to whole-brain, spatially distributed activity patterns, suggesting that estimated numerosity WM content is maintained in the LPFC (Nieder et al., 2002;Miller, 2003, 2004;Tudusciuc and Nieder, 2009;Spitzer et al., 2014a;Nieder, 2016).
It should be noted that we used temporally distributed tactile numerosity stimuli as the WM memoranda, namely the numerosity, was presented as a sequence of pulses. Evidence exists for potential differences in perceptual processing of spatially and temporally distributed numerosities, where spatially distributed stimuli appear to be processed in parietal regions while temporarily distributed stimuli do not (Cavdaroglu and Knops, 2019). In line with the finding of Cavdaroglu and Knops (2019), we used temporally distributed stimuli and did not find evidence of WM representations in posterior regions in our full brain FWE-corrected analysis. However, a small cluster (k ϭ 164) extending to right IPS was observed to represent remembered numerosity content at an uncorrected threshold of p Ͻ 0.001. While our results agree with numerosity WM findings in NHPs that suggest frontal rather than parietal coding for spatial numerosity stimuli during WM retention (for review, see Nieder, 2016), further investigation is needed to conclusively decide for the role of the IPS. The role of the IPS could be interpreted as specific to perceptual processing, and therefore was only revealed at a lower threshold in our analysis, while the PFC contains WM instead. Alternatively, a potentially different nature of the neuronal code (e.g. spatial distribution of a multivariate code) in the IPS might lead to the observed findings (Hebart and Baker, 2018). That is, it might be the temporarily distributed nature of the applied stimuli that drives the effects in the PFC, and the IPS would be more specialized for spatially distributed presentations as used by most previous studies. A future direct comparison of our results with spatial numerosity stimuli is necessary to test for differences determined by the stimulus types.
Moreover, while the literature relating to numerosity WM is limited, there is extensive work exploring the WM representation of abstract quantities more generally. Specifically, the frequency discrimination task has been systematically explored in a multitude of modalities with a wide range of methods (Romo et al., 1999;Lemus et al., 2009;Spitzer et al., 2010;Spitzer and Blankenburg, 2011;Fassihi et al., 2014;Vergara et al., 2016;Schmidt et al., 2017;von Lautz et al., 2017;Uluç et al., 2018;Wu et al., 2018). Numerosity and frequency share several traits, particularly that they are both abstract magnitudes that may be represented in a supramodal fashion (Nieder and Miller, 2003;Spitzer and Blankenburg, 2012;Nieder, 2016;Vergara et al., 2016). However, whether their underlying WM representations are maintained by a shared network has yet to be explored. The present study provides an initial step toward resolving this question by providing the first evidence that frequency and numerosity WM representations are maintained in overlapping brain regions. We identified numerosity-specific WM content in the right IFG, SMA, and left PMC, which is in agreement with results from frequency studies also using an fMRI MVPA approach in humans (Schmidt et al., 2017;Wu et al., 2018;Uluç et al., 2018). Unimodal and multimodal research in both NHPs and humans has identified frequency-specific content in the right LPFC and SMA, thereby suggesting that the WM representations are modality independent in nature (Romo et al., 1999;Hernández et al., 2002Hernández et al., , 2010Barak et al., 2010;Spitzer et al., 2010;Blankenburg, 2011, 2012;Vergara et al., 2016;Schmidt et al., 2017;Wu et al., 2018). However, the explicit relationship between frequency and numerosity still needs to be explored, particularly with respect to the underlying neural codes of numerosity and frequency representations (Nieder, 2017).
Additionally, we identified numerosity-specific content in the left PMC. Previous findings from frequency WM fMRI MVPA studies identified abstract quantity information in the PMC (Schmidt et al., 2017;Uluç et al., 2018;Wu et al., 2018) . Moreover, the dorsal PMC has been shown to represent abstract numerical rules, such as comparison and calculation (Gruber et al., 2001;Eger et al., 2003;Nieder, 2005). This is in line with the present task, which required the comparison of numerical quantities, suggesting representation of task-relevant, numerosity-specific information to be used in numerical comparison. All results are reported at p FWE-Cluster Ͻ 0.05 with a cluster-defining threshold of p Ͻ 0.001. Mean prediction accuracy over the delay period is reported. Areas were, where possible, identified using the SPM anatomy toolbox (Eickhoff et al., 2005).
In summary, the data at hand is in line with the suggestion of a domain general, abstract magnitude processing system. This abstract processing system can be identified by multivariate WM representations of tactile numerosity stimuli within the right PFC. Together with previous findings that found WM representations of tactile frequency (Spitzer et al., 2010(Spitzer et al., , 2014aSpitzer and Blankenburg, 2012;Schmidt et al., 2017;Wu et al., 2018), visual flicker frequency (Spitzer and Blankenburg, 2012;Spitzer et al., 2014a;Wu et al., 2018), auditory frequency Blankenburg, 2012, Uluç et al., 2018), and the reports of number coding (Nieder et al., 2002;Miller, 2003, 2004;Tudusciuc and Nieder, 2009;Nieder, 2016) in the PFC, the present study provides additional evidence suggesting that the PFC is capable of representing both analog quantities as well as parametric stimulus properties as frequencies. Thus, we provide preliminary evidence for a higher-level, modality-and format-independent, abstract quantitative WM system that resides within the PFC.