Tuned neural responses to haptic numerosity in the putamen

The ability to perceive the numerosity of items in the environment is critical for behavior of species across the evolutionary tree. Though the focus of studies of numerosity perception lays on the parietal and frontal cortices, the ability to perceive numerosity by a range of species suggests that subcortical nuclei may be implicated in the process. Recently, we have uncovered tuned neural responses to haptic numerosity in the human cortex. Here, we questioned whether subcortical nuclei are also engaged in perception of haptic numerosity. To that end, we utilized a task of haptic numerosity exploration, together with population receptive field model of numerosity selective responses measured at ultra-high field MRI (7T). We found tuned neural responses to haptic numerosity in the bilateral putamen. Similar to the cortex, the population receptive fields tuning width increased with numerosity. The tuned responses to numerosity in the putamen extend its role in cognition and propose that the motor-sensory loops of the putamen and basal ganglia might take an active part in numerosity perception and preparation for future action.

Many studies propose that the neurobiological foundation of numerosity is in a parietal-frontal cortical network ( Nieder and Dehaene, 2009 ). Electrophysiology studies in primates and crows revealed neurons that are tuned to numerosity in the parietal and frontal cortices of primates ( Nieder et al., 2002 ;Nieder and Miller, 2004 ), and in the nidopallium caudolaterale of crows, suggested to be the equivalent of the prefrontal cortex of primates Nieder, 2020 , 2016 ). In humans, similar tuned responses to numerosity were found in the intraparietal sulci ( Harvey and Dumoulin, 2017 ;Piazza et al., 2004 ), supporting the premise that numerosity is encoded by an approximate system ( Nieder and Dehaene, 2009 ). Though most neuroscientific studies of numerosity focus on the parietal and frontal cortices, the ability to perceive numerosity by a range of species suggests that some subcortical nuclei may be included in the process ( Nieder, 2020 ). Such an involvement is indicated by a recent psychophysical study that measured performance differences in monocular vs. dichoptic processing ( Collins et al., 2017 ). In addition, discriminating numerosity also elicits very rapid eye saccades, further providing behavioral evidence for the notion that numerosity is perceived by a primitive, automated system ( Castaldi et al., 2020 ). Collins and colleagues provided evidence for a ratio-dependent number processing system in the monocular portion of the visual system. These results, however, do not provide conclusive evidence and are limited in pointing at which subcortical areas take part in this processing system, especially given the notion that some extra-striate visual areas also process monocular information ( Rokers et al., 2009 ).
Here, we hypothesized that numerosity perception is not restricted to the cortex and that subcortical nuclei may partake in the active perception of numerosity. We have recently uncovered tuned neural responses to haptic numerosity in the human cortex ( Hofstetter et al., 2021 ). Build- In the equal individual size condition, each sphere was the same size irrespective of the number of spheres. In the equal total volume condition, the total volume of all spheres was kept constant. b. One to seven spheres were presented in ascending and descending order, interleaved by twenty spheres which served as baseline. This cycle was repeated twice in each functional run.
ing on these results, we utilized a task of haptic numerosity exploration, together with population receptive field (pRF) model of numerosity selective responses ( Harvey et al., 2013 ) measured at ultra high-field MRI (7T).

Participants
Seven volunteers participated in the study (1 female, 1 left handed, mean age 33, age range 26-46). All participants had normal haptic perception. All experimental procedures were cleared by the ethics committee of University Medical Centre Utrecht. Participants gave informed consent.

Haptic numerosity stimuli
We placed varying numbers of plastic spheres, suspended from plastic wires (1 to 7, with a baseline of 20 spheres; Fig. 1 ) in the right hand of participants. The spheres were all either the same size (equal size condition, where each individual sphere was 21 mm 3 ), or the same overall volume (equal volume condition, where the total sum of all spheres was 420 mm 3 ).
The haptic numerosity experiment included exploration periods of 3 s. Another set was placed in the hand after 4.5 s, i.e. the inter-stimulus interval. Participants laid in the scanner with their eyes closed and were asked to explore the spheres, but no numerosity judgments were required. The spheres were placed in the hands of the participant by an experimenter in the scanning room. The experimenter received auditory cues (through headphones the participant could not hear) indicating which numerosities to present at each point and when to place or remove the spheres. Numerosities 1 to 7 were presented in ascending order, followed by three presentation of 20 spheres. Then, numerosities 7 to 1 were presented in descending order, followed by a similar threetime presentation of 20 spheres ( Fig. 1 b ). This pattern was repeated 2 times during a functional run.
Previously we have shown that the perception of haptic numerosity in this type of experimental design follows the canonical responses of numerosity perception: quick an error-free judgment for numerosities in the subitizing range (up to four) and a significant increase in reaction time and judgment error for the higher numerosities. Moreover, measurements of fingers movements during this active task did not reveal any differences in the amount of unique motions, total amount of motion or the duration of motion with numerosity ( Hofstetter et al., 2021 ).
The numerosity experiments included 7-8 repeated runs in each scanning session. Each stimulus configuration (equal individual sphere size and equal total volume of spheres) was acquired in 1 or 2 scanning sessions on different days, so the overall each condition included 8 to 16 repetitions.

Functional data preprocessing
For each of the BOLD runs acquired per subject (across all conditions and sessions), the following preprocessing was performed: First, a reference volume and its skull-stripped version were generated using a custom methodology of fMRIPrep. A deformation field to correct for susceptibility distortions was estimated based on fMRIPrep's fieldmap-less approach. The deformation field is that resulting from coregistering the BOLD reference to the same-subject T1w-reference with its intensity inverted ( Huntenburg, 2014 ;Wang et al., 2017 ). Registration is performed with antsRegistration (ANTs 2.2.0), and the process regularized by constraining deformation to be nonzero only along the phase-encoding direction, and modulated with an average fieldmap template ( Treiber et al., 2016 ). Based on the estimated susceptibility distortion, a corrected EPI (echo-planar imaging) reference was calculated for a more accurate co-registration with the anatomical reference. The BOLD reference was then co-registered to the T1w reference using bbregister (FreeSurfer) which implements boundary-based registration ( Greve and Fischl, 2009 ). Co-registration was configured with six degrees of freedom. Head-motion parameters with respect to the BOLD reference (transformation matrices, and six corresponding rotation and translation parameters) are estimated before any spatiotemporal filtering using mcflirt (FSL 5.0.9, ( Jenkinson et al., 2002 )). The BOLD timeseries were resampled onto the following surfaces (FreeSurfer reconstruction nomenclature): fsnative, fsaverage. The BOLD time-series were The HRF parameters are estimated during the pRF fitting procedure. Different HRFs are found when fitting the model to the subcortex nuclei (blue) relative to the cortex (red). b. An example of pRF model predictions calculated using HRF parameters estimated in the subcortex nuclei (blue) and the cortex (red) and using the HRF default parameters (two-gamma). The measured time-series was taken from a voxel in the putamen. resampled onto their original, native space by applying a single, composite transform to correct for head-motion and susceptibility distortions. These resampled BOLD time-series will be referred to as preprocessed BOLD in original space, or just preprocessed BOLD. The BOLD timeseries were resampled into standard space, generating a preprocessed BOLD run in MNI152NLin2009cAsym space. First, a reference volume and its skull-stripped version were generated using a custom methodology of fMRIPrep. All resamplings can be performed with a single interpolation step by composing all the pertinent transformations (i.e. headmotion transform matrices, susceptibility distortion correction when available, and co-registrations to anatomical and output spaces). Gridded (volumetric) resamplings were performed using antsApplyTransforms (ANTs), configured with Lanczos interpolation to minimize the smoothing effects of other kernels ( Lanczos, 1964 ). Non-gridded (surface) resamplings were performed using mri_vol2surf (FreeSurfer).
The rest of the preprocessing was done in vistasoft ( https://github. com/vistalab/vistasoft/wiki ). The anatomical and functional data of each subject in MNI space was used (MNI152NLin2009cAsym). For each subject, the first 8 timeframes were discarded from the functional scans. To increase signal strength, data was also averaged across subjects. To that end, the data of each subject was first averaged across runs based on stimuli condition. Then, the averaged data of each subject was transformed into the anatomical template of MNI152NLin2009cAsym at a resolution of 1mm. No spatial or temporal smoothing was applied to the functional data. ROIs of the subcortical nuclei (caudate, putamen, globus pallidum, thalamus) and the ventricles were created using the MNI template (tpl-MNI152NLin2009cSym_res-1_atlas-CerebrA_dseg).

population receptive field modeling of numerosity
Numerosity responses in the subcortex were estimated using pRF modeling ( Dumoulin and Wandell, 2008 ;Harvey et al., 2013 ;Harvey and Dumoulin, 2017 ). Briefly, the pRF model describes the averaged tuning of the underlying neural populations using logarithmic Gaussian functions. These Gaussian functions are characterized by preferred numerosity (mean of the Gaussian) and tuning width (standard deviation of the Gaussian).
At each voxel, the pRF model is estimated based on the fMRI data and the time course of presented numerosities. For each candidate preferred numerosity and tuning width, a predicted neural response time course is calculated as the amplitude of the candidate neural response function at each time point's presented numerosity. Each candidate predicted neural response time course is then convolved with the haemodynamic response function (HRF) to create a candidate predicted fMRI time course. The chosen pRF parameters for each voxel are those whose predicted fMRI time course is best correlated with the voxel's measured fMRI time course. Next, we estimated the HRF parameters by comparing predicted time-series and actual time-series and chose the HRF parameters that minimized the difference between predicted and actual time series ( Harvey and Dumoulin, 2011 ). Last, the HRF parameters were used to re-fit the pRF.
The haemodynamic response function that was found to best predict the measured timeseries in the subcortical nuclei was faster than the canonical HRF typically used ( Glover, 1999 ) and from the HRF we have previously found to fit neural responses in the cortex ( Hofstetter et al., 2021 ) ( Fig. 2 ). This HRF pattern corresponds with former results of subcortical HRF at 7T MRI ( Lewis et al., 2018 ).
The pRF fitting procedure allows preferred numerosity estimates outside the range of the presented stimuli, ensuring estimates within the stimulus range are not just the best of a limited set. We excluded from analysis any recording sites where the preferred numerosity was outside our presented range (one to seven spheres).
To convert the variance explained (R 2 ) measures to probabilities of observing these model fits by chance, we built a null distribution using the same model-fitting procedure on 19,140 voxels taken from the ventricles from the data averaged across subjects and conditions. For each variance explained we calculated the proportion of these voxels with model fits exceeding that value. Subsequent analyses included only voxels in the subcortical nuclei where the variance explained of the pRF model was higher than 60% (equivalent to a probability of 0.0043 of observing this goodness of fit by chance). Similar analysis was done at the single subject level. Here, the null distribution was calculated using 133,980 voxels of the ventricles, gathered from data of all participants. The chance of the variance explained being higher than 60% at the single participant is p = 0.0005.

Cross validation
We distinguished between neural responses selective to numerosity and to the presence of stimuli in the hand and their subsequent exploration. In order to define recording sites where neural activity may be better explained by the motor component of the task, rather than tuning to stimulus numerosity, we split the data averaged across participants into two. These split data included data from both haptic conditions (equal total volume of spheres and equal individual spheres size) that were divided based on even-odd runs. We then fitted a general linear model (GLM) on each split data where the timing of stimulus presentation served as a predictor. We also fitted the numerosity pRF model on each half-split of data. Then, the goodness of fit of the pRF prediction from one half-split was evaluated on the time course of the other half, yielding new values of goodness of fit. Voxels where the variance explained by the GLM model was higher than the cross-validated pRF

Fig. 3. Neural responses to haptic numerosity exploration in the putamen. a.
Responses elicited by the haptic exploration of the numerosity stimuli and captured by the numerosity pRF model (variance explained). The variance explained by the numerosity model (R 2 ) is derived across the two stimulus conditions (see Fig. 1 ). b. Example of the fMRI time-series that illustrates similar responses throughout the task, irrespective of numerosity. The numerosities are indicated at the top of the graph. Due to the relatively slow stimulus presentation, the fMRI time-series reveals responses to each individual stimulus presentation. This pattern of responses was captured by a general linear model where the presence of the spheres serves as a predictor ( "on-off model ", solid line). c. Goodness of fit of numerosity-selective neural models after removal of points where the stimulus presentation better predicted the numerosity selective response predictions (see methods). Here, for illustration purposes, the variance explained is shown whithin the putamen. d,f. Two samples of fMRI time courses from sites in the putamen. The two sites show different responses to different haptic numerosities: (d) the largest response amplitude occurs at low numerosities; (f) the largest response occurs at higher numerosities. The numerosity pRF model's predictions (solid lines) capture these different responses as different numerosity selectivity parameters. e,g. The pRF model summarizes the neural responses using a logarithmic Gaussian tuning function with two parameters: preferred numerosity and tuning width, defined by the full width at half maximum (FWHM). numerosity models were excluded from further analysis. Similar analysis was performed at the single participant level. Here, data was split based on the even-odd runs taken from both haptic conditions. We excluded recording sites where the variance explained was lower than 60%, numerosity preference was outside stimuli range, or the fMRI time course was better explained by a response to stimuli presentation (and not to its numerosity, see Fig. 3 b ). We also fitted monotonic GLM to these data points that surpassed these criteria and compared the variance explained by the monotonic model to the pRF numerosity models.

Analysis of change in receptive fields size
We computed the tuning width as a function of preferred numerosity. The recording points were binned within every 0.25 increase in preferred numerosity. The mean of the tuning width and standard error of the mean within each bin was calculated.

Tuned responses to haptic numerosity in the putamen
Due to low SNR in the subcortex we averaged the functional data across participants. We then fitted a population receptive field (pRF) model of numerosity selective responses with two parameters: preferred numerosity and tuning width. The model was fitted to the subcortical nuclei (including the caudate, putamen, globus pallidum and thalamus; see supplementary Fig. 1a ). We found good fits of the pRF model in the putamen ( Fig. 3 a ). Next, we differentiated between neural responses selective to numerosity and to the exploration of the spheres, irrespective of the spheres numersoity (see methods ; Fig. 3 b ). We found that the pRF model explained a large amount of the signal variance in the putamen (mean variance explained = 73%). The pRF model revealed neural populations tuned to haptic numerosity in the putamen ( Fig. 3 c-g ). Last, we tested whether monotonic models may explain these neural responses better than the numerosity pRF models. However, the increasing or decreasing monotonic models showed a lower variance explained than the numerosity pRF models (averaged variance explained of increased and decreased monotonic models were 0.28 and 0.47, respectively).

Properties of haptic numerosity preferences in the putamen
Tuned responses to numerosity were found in the putamen. At the group level, numerosity preferences increased from anterior to posterior and medial to lateral in both nuclei ( Fig. 4 a ). This organization pattern was not influenced by one participant (see supplementary Fig. 2 ). Tuned responses to numerosity were also found at the single participant level, though to a lesser extent ( Fig. 4 b ). However, a clear organization of change in preferred numerosity across the putamen was only found using data averaged across participants.
Numerosity preferences that were obtained at each stimuli condition (i.e., equal total volume of spheres and equal individual sphere size) were well correlated (Pearson correlation, r = 0.42, p < 0.0001; Fig. 4 c  and supplementary Fig. 1b ). This correlation demonstrates that overall, numerosity preferences were similar across stimulus conditions and do not depend on the sphere's configuration.
We also found a link between tuning width and preferred numerosity. Similar to the tuned responses for haptic and visual numerosity that were found in the cortex ( Harvey et al., 2013 ;Harvey and Dumoulin, 2017 ;Hofstetter et al., 2021 ), the size of the numerosityselective receptive fields in the putamen increased with numerosity ( Fig. 4 d ), suggesting the receptive fields become less specific with increase in numerosity. Fig. 4. Haptic-driven numerosity preferences in the putamen. a. Numerosity preferences for data averaged across participants from both haptic stimulus conditions (variance explained > 60%). Colors represent preferred numerosity. b . Numerosity preferences for data averaged across both haptic stimulus conditions from two example participants. c . Pearson correlation analyses between preferred numerosity in the putamen of the two haptic stimulus configurations (i.e., equal total volume of spheres and equal individual sphere size). d . Increase in tuning widths plotted as a function of the preferred numerosity. Error bars show the standard error of the mean for each data point. Z, Y represents position in MNI space. For demonstration purposes, the results are shown whithin the putamen.

Summary of the main results
We found tuned responses to numerosity in the bilateral putamen when participants haptically explored spheres in their right hand. The size of the population receptive fields tuned to numerosity increased with increasing numerosity, similar to the pattern found in the cortex ( Harvey et al., 2013 ;Harvey and Dumoulin, 2017 ;Hofstetter et al., 2021 ) and in line with decreasing behavioral accuracy ( Hofstetter et al., 2021 ;Plaisier et al., 2009 ). The tuned responses to numerosity in the putamen correlated between stimuli conditions, indicating that the neu-ral responses to numerosity do not depend on the sensory configurations of the grasped objects. These results suggest that the putamen is involved in numerosity perception during active haptic exploration.

The putamen is involved in numerosity perception across modalities
Single case studies showed the involvement of the basal ganglia in numerical processing ( Benke et al., 2003 ;Corbett et al., 1986 ;Dehaene and Cohen, 1997 ;Delazer et al., 2004 ;Hittmair-Delazer et al., 1994 ) but the exact role of the different nuclei in the varied types of arithmetic calculations and mental processes is not clear. In healthy adults, imaging studies reported on activity in the putamen when participants learned to read Roman numerals ( Masataka et al., 2007 ), to compare between negative numbers ( Gullick et al., 2012 ), and through training on novel arithmetic problems ( Ischebeck et al., 2007 ). The activity in the putamen and caudate nucleus during tasks of negative numbers comparison and training in arithmetic calculations were suggested to be due to the striatum's role in rule and implicit learning ( Gullick et al., 2012 ;Ischebeck et al., 2007 ).
Only a few studies hinted at the involvement of the subcortical nuclei in approximate numerosity perception. Using a behavioral monocular/dichoptic paradigm, Collins and colleagues suggested that the subcortex takes part in numerosity discrimination ( Collins et al., 2017 ). This study, however, could not indicate which subcortical areas participate in the suggested monocular number processing system. Furthermore, behavior alone is not conclusive as some studies suggest that monocular information reaches extra-striate cortex ( Rokers et al., 2009 ). An imaging study of visual adaptation, where participants viewed the same approximate quantity repeatedly, found number-related decreasing activity in a network of brain regions in the cortex, cerebellum and the putamen ( Piazza et al., 2007 ). The putamen was also found among several cortical and cerebellar regions showing significant activation in a number comparison task relative to rest ( Göbel et al., 2004 ). These studies used visual stimuli for numerical comparisons or presentations. Our study, to the best of our knowledge, is the first to extend the question of numerosity perception in subcortical nuclei to the haptic modality and show that these subcortical nuclei display similar numerosity tuning properties as cortical regions ( Hofstetter et al., 2021 ). Our results propose that the putamen is involved in numerosity perception using haptic exploration, even when no numerical judgment is required.

Possible role of the putamen in the haptic numerosity system
As part of the basal ganglia-thalamocortical circuits, the putamen plays an important role in motor planning and execution ( Alexander et al., 1986 ). Dispersed input from similar cortical points is received in the putamen, where new sensory-motor associations may be formed before an output is sent to the globus pallidum ( Flaherty and Graybiel, 1994 ). Tracing studies in monkeys revealed that the putamen receives many cortical afferents ( Flaherty and Graybiel, 1994 ), amongst them are those from caudal frontal and rostral parietal lobes, including area 5 ( Alexander et al., 1986 ;Hedreen and Delong, 1991 ). Electrophysiological recording in monkey showed that neurons in this area respond to the number of self-generated movements ( Sawamura et al., 2010( Sawamura et al., , 2002. Similarly, direct connection to the rostral putamen were found from the medial temporal lobe (MTL), including the entorhinal cortex ( Hedreen and Delong, 1991 ). Tuned responses to symbolic and non-symbolic numerosities were recorded in the MTL in neurosurgical patients performing numerical calculations ( Kutter et al., 2018 ). The direct connectivity between the putamen to cortical areas involved in somatosensory and numerosity perception, as well as motor output, proposes that the uncovered tuned responses for numerosity in the putamen might be part of a sensory-motor system that links numerosity sensory perception with action ( Anobile et al., 2021 ). Such a link was found using behavioral adaptation techniques where finger tapping frequency affected perceived visual numerosity ( Anobile et al., 2016 ;Maldonado Moscoso et al., 2020 ). It should be noted, however, that fingers motion during the haptic exploration task that was used in the current study was similar across the explored numerosities ( Hofstetter et al., 2021 ). Therefore, the tuned responses to numerosity found in the putamen are not suggested to reflect the numerosity of specific self-generated movements. We hypothesize that tuned responses to numerosity in the putamen may be part of the putamen's role within the sensory-motor network, where numerical information, such as the number of items held by the hand, might be relevant for future movements and actions.

Study limitations
The current study aimed at investigating the involvement of subcortical nuclei in numerosity perception. However, the study of the basal ganglia and small nuclei of the brainstem with functional MRI at ultra-highfield bears more technical limitations than imaging BOLD responses at cortical regions. Lower BOLD sensitivity in the basal ganglia is caused by higher noise levels, reduced T2 * relaxation values, and lower signal levels due to distance from the receiver head coil ( Aquino et al., 2009 ;de Hollander et al., 2017 ;van der Zwaag et al., 2016 ). As a result, our ability to detect reliable tuned reposes to numerosity at the single subject level were limited and were found in only a few participants.
We have also analyzed a previous dataset of visual numerosity perception ( Hofstetter et al., 2021 ) but did not find evidence in support of tuned neural responses as the variance explained by the model was much lower and did not survive the 60% variance explained threshold. This could mean that the neural populations in the putamen are specific for haptic numerosity only. However, as mentioned above, a few studies mentioned some responses in the putamen during tasks of numerical perception using visual stimuli ( Göbel et al., 2004 ;Piazza et al., 2007 ). Furthermore, our failure to find neural populations tuned to visual numerosity does not mean they do not exist. Future studies would benefit both from more data acquisition and from protocols that are tailored for imaging the subcortical nuclei.
Due to the limited SNR in the subcortex we have averaged data across participants. Averaging across participants might bias the peak of the estimated receptive fields towards the center of the stimulated range and cause bigger tuning widths due to a larger variability of the underlying neural responses. However, the results of the current study demonstrate that even in an averaged data we were not only able to detect tuned responses to the entire range of the presented stimuli, but the tuned responses showed both a spatial organization across the putamen as well as a relation between tuning width size and numerosity preferences.

Conclusions
Many species, including humans, share an intuitive capacity for numerosity perception. Correct numerosity estimation is critical when formulating a correct behavioral response ( Nieder, 2020 ). In the current study we have found tuned responses to haptic numerosity in the bilateral putamen. These results extend the role of the putamen in cognition and propose that the motor-sensory loops of the putamen and basal ganglia might take an active part in numerosity perception and preparation for future action.

Declaration of Competing Interest
There is no conflict of interest to declare.

Data and code availability
The code that supports the findings of this study is available from the corresponding author upon reasonable request.

Funding
This work was supported by the Ammodo KNAW Award to S.D. and the Netherlands Organization for Scientific Research grant (NWO-VICI 016.Vici.185.050) to S.D.

Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.neuroimage.2021.118178 .