A touch of hierarchy: population receptive fields reveal fingertip integration in Brodmann areas in human primary somatosensory cortex

Several neuroimaging studies have shown the somatotopy of body part representations in primary somatosensory cortex (S1), but the functional hierarchy of distinct subregions in human S1 has not been adequately addressed. The current study investigates the functional hierarchy of cyto-architectonically distinct regions, Brodmann areas BA3, BA1, and BA2, in human S1. During functional MRI experiments, we presented participants with vibrotactile stimulation of the fingertips at three different vibration frequencies. Using population Receptive Field (pRF) modeling of the fMRI BOLD activity, we identified the hand region in S1 and the somatotopy of the fingertips. For each voxel, the pRF center indicates the finger that most effectively drives the BOLD signal, and the pRF size measures the spatial somatic pooling of fingertips. We find a systematic relationship of pRF sizes from lower-order areas to higher-order areas. Specifically, we found that pRF sizes are smallest in BA3, increase slightly towards BA1, and are largest in BA2, paralleling the increase in visual receptive field size as one ascends the visual hierarchy. Additionally, we find that the time-to-peak of the hemodynamic response in BA3 is roughly 0.5 s earlier compared to BA1 and BA2, further supporting the notion of a functional hierarchy of subregions in S1. These results were obtained during stimulation of different mechanoreceptors, suggesting that different afferent fibers leading up to S1 feed into the same cortical hierarchy.


Introduction
Touch is an important source of information about our direct surroundings. We use touch information to explore objects and surfaces and touch plays a major part in haptic processes such as tool use. The loss of adequate touch signal processing, e.g. due to stroke, frequently leads to severe impairments affecting many facets of everyday life. Hence, understanding somatosensory processes in the human brain following cutaneous touch signals is relevant to many scientific areas ranging from fundamental neuroscience to the deciphering of neurological disorders of the somatosensory system. Imaging studies in humans have mostly addressed the somatotopic organization of the hand and fingers (Maldjian et al. 1999;Kurth et al. 2000;Hlustík 2001;Blankenburg et al. 2003;Nelson and Chen 2008;Schweizer et al. 2008;Sanchez-Panchuelo et al. 2010;Ann Stringer et al. 2014;Martuzzi et al. 2014;Choi et al. 2016;Kikkert et al. 2016;Kolasinski et al. 2016;Sanchez Panchuelo et al. 2018;Rocha et al. 2019;Puckett et al. 2020), and the whole body (Akselrod et al. 2017;Tal et al. 2017). However, other functional characteristics of human S1 have not received equal attention. Specifically, the processing hierarchy of cyto-architectonically distinct regions in human S1, i.e. Brodmann areas BA3a/b, BA1, and BA2, (Brodmann 1909;Geyer et al. 1999), has been investigated structurally in humans (Sánchez-Panchuelo et al. 2014;Wagstyl et al. 2015), but not from a functional perspective. In the current study, we investigate the functional hierarchy in human S1 1 3 by estimating the integration of somatic information in different Brodmann areas.
When cortical information is processed at different hierarchical levels, information from multiple lower-level sources is integrated at the higher-order level. As a result, regions of higher hierarchical order contain neurons that exhibit larger or more complex receptive fields, meaning that neurons are responsive to more input or specific combinations of input. Functional hierarchy among separate S1 regions in humans can, therefore, potentially be revealed through a form of spatial somatosensory information integration (Hubel and Wiesel 1968;Duffy and Burchfiel 1971;Van Essen and Maunsell 1983). Previous animal studies have reported that BA3b is the primary target of thalamic output from the ventrolateral and ventroposterior nucleus (Jones and Powell 1970;Chung et al. 1986;Miller et al. 2001), which then projects onwards to BA1 and BA2 (Friedman 1983;Felleman and Van Essen 1991;Kaas 1993;Iwamura 1998). As a result, neuronal receptive fields, as reported in animal studies, are smallest in BA3b and increase in size in BA1, BA2 and beyond (Armstrong-James 1975;Hyvärinen and Poranan 1978;Sur et al. 1980;DiCarlo et al. 1998). In humans, receptive field properties of individual neurons cannot easily be assessed in healthy volunteers under normal circumstances. However, average receptive field properties of small neuronal populations (e.g. neurons inside a single MRI-voxel) can be estimated using a Gaussian population Receptive Field (pRF) model. PRF modeling was originally developed for vision (Dumoulin and Wandell 2008), where it has exposed hierarchical processing characteristics as well as other traits of the human visual system (Harvey and Dumoulin 2011;Haak et al. 2012;Dumoulin et al. 2014;Klein et al. 2014;Wandell and Winawer 2015;Merkel et al. 2018;Welbourne et al. 2018). Furthermore, two recent functional MRI (fMRI) studies have shown that pRF modeling can also be used to describe the average receptive field properties of small neuronal populations in human S1 (Schellekens et al. 2018;Puckett et al. 2020). Even though some studies find evidence consistent with hierarchical organization of somatosensory processing in humans (Bodegård et al. 2001;Van Boven et al. 2005;Dijkerman and de Haan 2007;Kim et al. 2015;Whitehead et al. 2019), the extent of spatial integration across different Brodmann areas in human S1 is presently not well defined.
The current objective is to estimate pRF properties across Brodmann areas, following vibrotactile stimulation of the fingertips. Vibrotactile stimulation can be signaled by two distinct cutaneous mechanoreceptors: Meissner corpuscles and Pacinian corpuscles, depending on the frequency of vibration (Mountcastle et al. 1972;Bolanowski et al. 1988;Pasterkamp 1999). Meissner corpuscles typically show a peak activity for flutter frequencies (i.e. between 10 and 50 Hz), while Pacinian corpuscles respond to higher frequencies with a preference around 250 Hz (Rowe 2002). Furthermore, previous studies showed that Meissner and Pacinian corpuscles signal somatosensory information through different pathways, i.e. Rapid-Adapting (RA) and Pacinian pathways (Vallbo and Johansson 1984;Gescheider et al. 2004;Harvey et al. 2013;Saal et al. 2015), which reportedly project to different regions of the thalamus (Herron and Dykes 1986;Kaas 1993). Additionally, Pacinian pathways may have more connections to BA1 than BA3b (Paul et al. 1972;Hyvärinen and Poranan 1978;Iwamura et al. 1993). Hence, the hierarchical order of somatosensory processing among Brodmann areas in S1 may be frequency-dependent or at least influenced by the supplied frequency of vibration. To investigate hierarchical differences caused by stimulated mechanoreceptor type, we supplied a vibrational stimulus to the fingertips at three different frequencies: 30 Hz, 110 Hz, and 190 Hz. A perfect isolation of stimulated mechanoreceptor type is not realistic and multiple pathways likely contribute to the observed cortical signal with increasing contributions of Pacinian pathways for higher stimulation frequencies (Choi et al. 2016;Kuroki et al. 2017). Thus, differences in initial cortical projection site between RA an Pacinian pathways could be detected through changes in pRF size for different vibrotactile stimulation frequencies.
In the present study, we scrutinize the hierarchical organization of S1 by measuring the properties of tactile pRFs in BA3b (from here on referred to as BA3), BA1, and BA2. The five fingers of the right hand were vibrotactually stimulated at three different frequencies, 30 Hz, 110 Hz, and 190 Hz, while Blood-Oxygen-Level-Dependent (BOLD) activity in S1 was measured with 7 T fMRI. PRF modeling allows us to infer the somatotopic tuning of neuronal populations in each of the three Brodmann areas. We expect an increase in pRF size, the specificity of the somatotopic tuning, along the somatosensory processing pathway. Such a finding would indicate increasing spatial integration and be in accordance with sequential information processing and increasing processing complexity from BA3 to BA1, and finally BA2. The hierarchical order across Brodmann areas is further investigated by examining the temporal dynamics of the hemodynamic response function (HRF). Finally, the effect of mechanoreceptor pathway on cortical pRF size is presently unknown. Through pRF size estimations in different Brodmann areas under different vibrotactile frequency conditions, we investigate putative differences in cortical hierarchical projections related to different mechanoreceptor types.

Participants
Eight healthy volunteers (age range 23-31 years old, 4 female) participated in the study. All participants gave written informed consent before entering the study. The protocol was approved by the local medical ethics committee of the University Medical Center Utrecht, Netherlands, in accordance with the Declaration of Helsinki (2013).

Apparatus
The vibrotactile stimulus was delivered using MR-compatible piezoelectric stimulators with a triangular shaped tip and a contact area of approximately 1 mm 2 (http:// dance rdesi gn. co. uk/). The stimulation was controlled via a custom-written MATLAB (www. mathw orks. com) script. Analog stimulus signals were transferred to the stimulators using a NI-9264 digital-to-analog converter output module (National Instruments, Austin, TX, USA), which was connected to a conventional laptop and an amplifier.
We mounted five stimulators on a plexiglass plate using ordinary adhesive gum. The adhesive gum allowed for the repositioning of the five stimulators to match each participant's hand. The fingertips of the right hand were placed on the stimulators (digits did not touch each other). The hand and fingers were taped to the plexiglass plate with standard paper tape to prevent the fingers from accidentally disconnecting from the stimulators. The plexiglass plate rested on the participant's abdomen, while the right elbow was supported by towels. Using this setup, the subject could maintain a stationary position of the right arm/hand comfortably for the full length of the fMRI experiments. This minimized movement of the hands, which could affect the results. Moreover, subjects were explicitly instructed to keep both hands still during the experiments.

Procedure and stimuli
Each subject underwent 4 fMRI experiments: the first 3 were pRF experiments, conducted to estimate pRF properties (i.e. receptive field center, size, and amplitude). These 3 experiments differed only with respect to the frequency of vibration (30 Hz,110 Hz,and 190 Hz). The 4th fMRI experiment was conducted to estimate the hemodynamic response function (HRF) within each individual subject's S1. During the 3 pRF experiments, each fingertip was stimulated 8 times in a pseudo-randomized order. Only one fingertip was stimulated at a time, and a single stimulation lasted for 4 s. An intermittent stimulation paradigm was chosen to minimize adaptation processes and, therefore, maximize the observed BOLD response: during the 4 s stimulation period, a 400 ms on period was alternated with a 100 ms off period. After the 4 s stimulation period, a 10 s rest period ensued except for 8 randomly selected stimulation periods when the ensuing rest period was lengthened to 14.4 s. Our analysis did not require a complete return to baseline, but rather allowed for the response to one stimulus to persist into the onset of the next. In total, a single pRF experiment took 595.2 s. During the HRF experiment, a brief vibrotactile stimulation of 500 ms at 30 Hz was applied to all 5 fingertips simultaneously. The brief 500 ms stimulation was delivered intermittently: 200 ms on/100 ms off/200 ms on. There were 32 500 ms events throughout the HRF experiment with variable inter-stimulus interval (ISI). The minimum ISI was 3.05 s, the maximum ISI was 23.97 s, and the median ISI was 7.98 s.The full HRF experiment took 320 s.

Scan protocol
Scanning was conducted at a 7 Tesla Philips Achieva scanner (Philips, Best, Netherlands), using a volume transmit and a 32-channel receive headcoil (Nova medical, MA, USA). A multi-slice gradient echo (GE) echo-planar imaging (EPI) sequence was used for functional image acquisition with the following specifications: TR/TE: 1600/27 ms, flip angle: 70°, SENSE factor: 3 in the anterior-posterior direction, field-of view (FOV) (ap,fh,rl): 209.4 × 41.6 × 165.0 mm at 1.6 × 1.6 × 1.6 mm voxel resolution, and interleaved slice acquisition. The FOV was placed on the superior part of the brain, covering the hand region of the postcentral gyrus. 372 volumes were acquired per pRF experiment and 200 volumes were acquired for the HRF experiment. Additionally, 10 volumes were acquired with a reversed phase encoding direction (i.e. posterior to anterior) for correction of geometrical distortions. Finally, a whole-brain T1-weighted volume was acquired with TR/TE: 7.00/3.05 ms, flip angle: 8°, FOV (ap,fh,rl): 250 × 200 × 190 mm at 0.78 × 0.78 × 0.8 mm voxel size, and a whole-brain proton density volume of equal dimensions.

Image processing
The T1-weighted anatomical volume was adjusted for proton density to correct for large scale intensity inhomogeneities (Van de Moortele et al. 2009). Afterwards white matter and pial brain surfaces were estimated using Freesurfer (https:// surfer. nmr. mgh. harva rd. edu/). These surfaces were also inflated and flattened using Freesurfer. The functional volumes were slice time corrected, realigned (i.e. corrected for head motion), corrected for geometrical distortions, and co-registered to the anatomical T1-weighted volume using AFNI. Transformation matrices for these steps were computed using the AFNI functions 3dvolreg, 3dQwarp, and 3dAllineate, respectively. The transformation matrices were combined and all spatial preprocessing transformations were applied within a single interpolation step using the AFNI function 3dNwarpApply to minimize smoothing caused by multiple interpolation steps and general interpolation errors. The functional volumes were mapped onto the estimated cortical surface reconstructions across the full depth of the 1 3 estimated gray matter using Freesurfer, creating a timeseries per surface vertex. The timeseries were high-pass filtered with a cut-off at 0.01 Hz and rescaled to percent signal change. Finally, regions of interest were drawn on the reconstructed cortical surface, based on the Brodmann area atlas supplied by Freesurfer (Fischl et al. 2008). Region BA3 corresponded with atlas areas BA3a and BA3b (covering the rostral wall of the postcentral gyrus). Region BA1 corresponded with atlas area BA1 (covering the crown of the postcentral gyrus). Finally, region BA2 (covering the caudal wall of the postcentral gyrus) was based on atlas area BA2, but manually limited posteriorly at the base of the postcentral sulcus.

pRF analysis
Each vertex' timeseries was fitted with a Gaussian receptive field model, which described the signal amplitude for any fingertip stimulation (1): where "x i " represents the stimulated fingertip and "N" is the list of fingertips ranging from 1 = thumb to 5 = little finger. The estimated pRF center, "x 0 ", describes the preferred fingertip per surface vertex and can be any real number (including fractioned numbers) between 0.5 and 5.5. A surface vertex is taken to prefer: the thumb when, 0.5 < "x 0 " < 1.5, index finger when, 1.5 < "x 0 " < 2.5, middle finger when, 2.5 < "x 0 " < 3.5, ring finger when, 3.5 < "x 0 " < 4.5, and the little finger when, 4.5 < "x 0 " < 5.5. The estimated pRF size, "σ", is the spread of the Gaussian in units of fingers: the larger the pRF size, the more the neuronal population responds to stimulated fingertips in addition to the preferred one. The receptive field model "g(x i )", then, is used to construct the effective task design (2): where "r(t)" is the effective task design, "s(x i ,t)" is the onset design matrix, which is a 2D binary matrix representing for each fingertip "x i " the stimulation onset and duration in scans "t". The multiplication of the onset design matrix "s(x i ,t)" and the Gaussian receptive field model "g(x i )" is summed over the fingertip dimension, resulting in the effective task design "r(t)". The effective task design is convolved with a hemodynamic response function (HRF), resulting in the predicted timeseries (3): (1) where "h(t)" is the HRF. Instead of assuming a canonical HRF, we convolved the estimated HRFs from the HRF experiment (averaged across subjects, see below) with the effective task design "r(t)". Therefore, we used an HRF that was specific for each Brodmann area. The predicted timeseries model "p(t)" was compared with the measured timeseries of each vertex (4): where y(t) is the measured vertex' timeseries, "p(t)" is the predicted timeseries, "β" is a scalar representing the signal amplitude and "c" is a constant. During the fitting procedure, optimal fits are calculated for the pRF center "x 0 " and size "σ" from Eq. (1) and "β" and "c" from Eq. (4) using the Levenberg-Marquardt (Markwardt 2009) least-square minimization algorithm (Fig. 1). Finally, goodness-of-fit F-statistics were calculated for each surface vertex model fit.

HRF analysis
For the HRF experiment, we estimated the hemodynamic response function of each vertex using a set of finite impulse response (FIR) functions (Lindquist et al. 2009). The timeseries were upsampled by a factor of four using a 3th degree B-spline interpolation, resulting in a time point every 400 ms. This matched the stimulus onset resolution, as stimulus onsets were locked to time samples every 400 ms. A set of finite impulses were constructed to cover the range of 14.4 s (i.e. 36 finite impulses), starting from the moment of stimulation. The amplitude in percent signal change at each time point was calculated using a multiple linear regression. An HRF per ROI was created by averaging the estimated HRFs of all vertices within the ROIs that showed a significant fit with respect to the HRF task design (false-discovery-rate corrected). Afterwards, the peak amplitude, time to peak (TTP) and full-width-at-half-maximum (FWHM) were extracted from the estimated HRF curves.

Statistical analyses
For the statistical analyses of all experiments, we included the surface vertices with a significant goodness-of-fit F-statistic derived from the pRF experiments (false-discoveryrate corrected) that fell in one of the three predefined ROIs. The percentage explained variance per vertex was calculated through the Pearson correlation coefficient of predicted timeseries and obtained timeseries squared. The presence of a somatotopy was assessed using the vertex coordinates of the flattened surfaces. Initially, the flattened surfaces were manually rotated so that the central sulcus was vertically aligned along the dorsoventral axis. A somatotopy is defined here as the linear relationship between dorsoventral (4) y(t) = • p(t) + c coordinates and pRF centers. Hence, the slope between dorsoventral coordinates and pRF centers reflects the presence of a somatotopy, given in pRF center per mm flattened surface, and was calculated using a linear regression per ROI, per vibrotactile frequency, and per subject. We used Student's t-test to test if slopes deviated significantly from zero. We used a 2-way univariate repeated measures ANOVA with the slopes as dependent variable and ROI and vibrotactile frequency as repeated measures factors (3 levels each) to test for differences in somatotopic structures per ROI or frequency of vibration. The pRF sizes were binned in five preferred finger representation bins, according to the pRF centers. Then, we applied a 3-way univariate repeated measures ANOVA to test for differences in pRF size across ROI, vibrotactile frequency, and preferred finger representation (with 3, 3, and 5 levels, respectively) with linear contrasts for each factor. The same 3-way univariate repeated measures ANOVA was performed on the estimated amplitude of the percent BOLD signal change (i.e. "β" from Eq. (4)). For the HRF experiment, differences in peak amplitude, TTP, and FWHM per ROI were also tested for using univariate repeated measures ANOVAs with only ROI as factor (3 levels). Additionally, we conducted the full pRF analysis using a canonical HRF for comparison purposes. We used a paired sample t test to compare the somatotopy slopes of the Brodmann area-specific HRFs with the canonical HRF.

Simulation analysis
Finally, we performed a simulation analysis to test for the influence of noise on estimated pRF parameters. We constructed model pRF timeseries on the basis of all possible combinations of pRF parameters (see pRF analysis), convolved with each of the 3 ROI-specific HRFs and added For visibility, only a part of the complete timeseries is shown. The onsets of the fingertip stimulation conditions are represented by the colored bars, see also hand icon. This particular vertex was acquired from subject 4, BA1, 190 Hz, and was fitted with a model with pRF center = 2.74 (between index and middle finger) and pRF size = 1.70 finger units random normally distributed noise to these model pRF timeseries. The added noise was equal in magnitude to the estimated noise from the original fMRI data set, which we estimated as the standard deviation of all included surface vertices' timeseries from all participants, after subtraction of the pRF model fit (i.e. residuals). Then, pRF parameters were estimated from the pRF model timeseries with added noise, and compared to the original noise-free pRF model parameters. The comparison was calculated as the percentage deviation of noisy parameters from noise-free parameters including the 95% confidence interval (95% CI). This procedure was iterated 100,000 times for each ROI-specific HRF, resulting in 300,000 simulated pRF timeseries.

S1 Somatotopy-spatial organization of pRFs
We used a Gaussian receptive field model to estimate the timeseries of the pRF experiments (Fig. 1B). The predicted timeseries explained on average 35% (s.d. = 11%) of variance of the recorded BOLD fMRI signal within the 3 predefined ROIs. On the basis of the estimated pRF centers we found the somatotopy of the five fingertips along the ventrolateral to mediodorsal axis of the postcentral gyrus in all 3 Brodmann areas (Fig. 2): BA3: t (7) = 13.10, p < 0.001, BA1: t (7) = 13.25, p < 0.001, BA2: t (7) = 8.51, p < 0.001. The somatotopy, characterized as the slope of cortical coordinates and pRF centers, differed significantly across the 3 Brodmann areas (F (2,14) = 15.26, p < 0.001). Particularly, the somatotopy was less clear in Brodmann area BA2 (post-hoc somatotopy slope t tests BA3-BA1: t (7) = 0.55, p = 0.589; BA3-BA2: t (7) = 5.04, p < 0.001; BA1-BA2: t (7) = 4.48, p = 0.001. This effect did not change when using a canonical HRF (t (8) = 0.71, p = 0.499), meaning that any observed somatotopy is not likely affected by the selected HRF. In BA2, there appears to be a cluster of pRF centers for the thumb and index finger and a second cluster for the middle, ring and little fingers (Fig. 2B). The frequency of vibration, however, did not influence the somatotopy slope (F (2,14) = 0.25, p = 0.782), although the projected somatotopy appeared less clear in several participants during the 30 Hz vibration condition compared to higher frequencies (Fig. 3). We, finally, did not observe an interaction effect between Brodmann areas and applied frequency of vibration on the somatotopy slope (F (4,28) = 0.85, p = 0.505), meaning that we did not find evidence for a somatotopy change in any Brodmann area for higher frequencies.

Amplitude of the BOLD signal
We found that the amplitude of the estimated percentage of BOLD signal change ("β", Eq. (4)) differed significantly across the 3 Brodmann areas (F (2,14) = 8.15, p = 0.004), where largest percent signal changes were measured in BA3 and gradually decreased towards BA2 (t (14) = − 4.03, p = 0.001, Fig. 5D). However, both preferred fingertip and vibrotactile frequency did not have a significant effect on the BOLD signal amplitudes (F (4,28) = 2.21, p = 0.094, and F (2,14) = 1.75, p = 0.208, respectively, Fig. 5E-F). Thus, the percent BOLD signal change differed per Brodmann area, but was not significantly affected by the preferred fingertip of included populations, or by the vibrotactile frequency at which fingertips were stimulated.
Since we found that both pRF size and signal amplitude differed across Brodmann areas (Fig. 5A, D), a methodological concern is that differences in signal-to-noise-ratio (SNR) between Brodmann areas could explain the differences in pRF estimates. Therefore, we performed a post-hoc simulation analysis to test for the influence of noise on pRF size and signal amplitude. The simulations indicated that the signal amplitude increased by approximately 11% (95% CI 10.6-11.6%) and the pRF size decreased by 1.5% (95% CI 1.3-1.7%) when noise-equal in magnitude to noise in the original fMRI data set-was added to the pRF model. This result indicates that noise may have influenced the estimates of signal amplitudes across Brodmann areas, which were measured to be in the order of 15% (BA3-BA1) to 25% (BA1-BA2), whereas it likely had a minor effect on the pRF size differences across Brodmann areas, which we found to be in the order of 20% (BA3-BA1) to 60% (BA1-BA2).

Fig. 2 Fingertip somatotopy.
A Single subject pRF centers following 190 Hz vibrotactile stimulation are presented on a pial surface and flattened surface (circle). The cortical coordinates along the dorsoventral axis plotted against the pRF centers are shown for all three Brodmann areas. For the pRF centers, 1 = thumb, 2 = index finger, 3 = middle finger, 4 = ring finger, 5 = little finger, which is also indicated by the colors in the scatterplot and the hand icon. B Group average of cortical coordinates along the dorsoventral axis plotted against the mean pRF center per fingertip 1 = thumb, 2 = index finger, 3 = middle finger, 4 = ring finger, 5 = little finger). Shaded area represents standard error of the mean across subjects. Different symbols represent different vibrational frequencies

General discussion
In the current study, we estimated pRFs in 3 subdivisions of human S1. The patterns of pRFs can be used to suggest a cortical hierarchy among these areas, if we operationalize the notion of hierarchy by the size of receptive field, specifically assuming that an area with smaller pRFs is earlier in the hierarchy. We fitted a pRF model to fMRI BOLD activity in S1, following vibrotactile stimulation of the fingertips. Additionally, we stimulated at 3 different frequencies of vibration to investigate changes in pRF size across S1 related to mechanoreceptor type and corresponding afferents. We found that pRF sizes increased from BA3 to BA1 and finally BA2, consistent with the notion of a cortical hierarchy in which spatial somatic information is pooled into larger and larger regions. This effect was observed under all vibrotactile frequency conditions. PRF sizes also increased with higher frequency of stimulation. These latter two results suggests that RA and Pacinian channels share a similar cortical hierarchy, but that somatic information from a relatively larger area of the hand is pooled in S1 neuronal populations during stimulation at higher frequencies. During all frequencies  The base of the central sulcus is shown by the white downward triangle, and the crown of the postcentral gyrus is indicated by the black upward triangle of vibrotactile stimulation we observed a somatotopy of fingertips, despite the somatotopy being less clear in BA2 compared to BA3 and BA1. No significant effect of frequency on somatotopy was observed, indicating that the whole of S1 responds to vibrotactile fingertip stimulation regardless of stimulation frequency. Finally, we found that pRF sizes gradually increased from thumb to little finger. Neuronal populations that preferentially code for the thumb responded least to stimulation of other digits, compared to neuronal populations coding for the little finger, which responded to stimulation of most other digits.

Cortical hierarchy S1
Cortical hierarchy was defined in this study through information integration, which increases when information progresses higher up the processing hierarchy. Information integration is associated with the widening of response profiles  Figure  shows the average pRF size across subjects for Brodmann areas A, fingertip representation B, and vibrotactile frequency C, as well as the corresponding estimated BOLD signal amplitude (D-F). Error bars denote the standard error of the mean across subjects of neuronal populations with respect to information coming from any number of possible sources. We estimated the widening of the response profiles of neuronal populations with a Gaussian shaped population receptive field model, where the spatial integration of somatosensory information is represented by the pRF size. We find that pRF sizes differ substantially between Brodmann areas, BA3, BA1, and BA2. Neuronal populations in BA3 have on average smallest pRF sizes, and the pRF sizes increase along the cortical processing hierarchy towards BA1 and are largest in BA2. PRF sizes in BA2 are approximately twice the size as the pRF sizes measured in BA3, of which a mere 1.5% can be explained by differences in SNR This result is likely analogous to the pRF size increase among cortical areas in visual cortex, where the primary visual cortex (V1) predominantly receives thalamic output and exhibits smaller receptive field sizes than visual cortical areas further up the hierarchy, as measured both at the single unit level (Felleman and Van Essen 1991) and the population level with fMRI (Dumoulin and Wandell 2008;Wandell and Winawer 2015), which likely reflects the average receptive field size of the underlying ensemble of neurons.
The hierarchical order of BA3, BA1 and BA2 is further supported by a shorter time-to-peak of the estimated HRF in BA3 compared to BA1 and BA2, which has also been observed in magnetoencephalography (MEG) studies (Inui et al. 2004;Suzuki et al. 2013). Thus, the order of cortical processing becomes apparent not merely through information integration, but also in the temporal domain. However, it is important to note that both feedforward and feedback neuronal processes contribute to the observed HRFs. Therefore, differences in temporal components of the HRF cannot solely be attributed to differences in sequential processing order. It is, for instance, possible that populations in BA1 and BA2 are not merely involved in somatosensory processing at a later point in time, but also for a slightly longer period of time, which would influence the observed HRF. Additionally, HRF latency can be affected by nonneural processes, such as the presence of draining veins (Lee et al. 1995). Nevertheless, the time-to-peak of the observed HRF in BA3 is roughly 0.5 s faster compared to the time-topeak of the HRF in BA1 and BA2. Assuming factors such as draining veins do not vary systematically between subareas in S1, this difference likely has a neuronal contribution. Our findings extend animal findings to humans, and are consistent with a cortical hierarchy in human S1, in which BA3 is the first cortical area to receive tactile information, which is then forwarded to BA1 and BA2.

Mechanoreceptive afferents
We applied three frequencies of vibrotactile stimulation to the fingertips to investigate the cortical hierarchy in human S1 as a result of different cutaneous mechanoreceptor afferents. The 30 Hz flutter frequency most likely activated Meissner corpuscles, whereas the higher frequencies would have resulted in increased contributions of Pacinian corpuscles (Bolanowski et al. 1988;Johnson 2001). Regardless of the stimulated mechanoreceptor, we observed somatotopic structures in all three included Brodmann areas. However, the somatotopy in BA2 was less clear than in the other two areas, which likely reflects less clear distinctions between cortical finger representations for areas higher up the cortical hierarchy, which has been reported in a previous animal study (Iwamura et al. 1983(Iwamura et al. , 1993Pons et al. 1985). We did not observe that the frequency of vibrotactile stimulation influenced the somatotopic structures of Brodmann areas, which may be in agreement with the notion of S1 neurons responding to multiple mechanoreceptor modalities (Pei et al. 2009;Abraira and Ginty 2013;Saal and Bensmaia 2014). However, previous optical imaging studies in monkeys have observed distinct columnar structures related to different types of mechanoreceptors in BA3. (Chen et al. 2001;Friedman et al. 2004). These frequency-dependent cortical columns are reportedly smaller than 400 µm in size. The spatial resolution used in this study was not sufficiently high to capture these differences in cortical projection for different mechanoreceptor afferents.
Our results show that pRF sizes increase with increasing frequency of vibrotactile stimulation. This effect was not found to differ across the three Brodmann areas and, therefore, we find no evidence to support the notion that different mechanoreceptors types project to S1 in different ways. The increase in pRF size for increased frequency Fig. 6 Hemodynamic response functions. Estimated hemodynamic response functions per Brodmann area. The areas denote one standard error of the mean across subjects could have been caused by several different processes. First, cutaneous mechanoreceptive units have receptive fields themselves, which could shape the feedforward information stream to S1. Mechanoreceptors in glabrous skin such as the Meissner corpuscle have relatively small receptive fields, whereas Pacinian corpuscles reportedly have receptive fields that extend beyond the range of one finger (Bell et al. 1994;Bolanowski and Pawson 2003). Second, neuronal activation thresholds could be dependent on vibrotactile frequency (Nelson et al. 2004;Simons et al. 2005;Ryun et al. 2017). Suprathreshold levels of activity for S1 neuronal populations could be attained during stimulation of cutaneous mechanoreceptors at high frequencies that would fall outside the neuronal populations' receptive fields during stimulation at lower frequencies. Third, the increase in the observed pRF size for higher frequencies of vibrotactile stimulation might be an extra-classical receptive field effect (Friston 2005;Schwabe et al. 2006). It has been suggested that vibrotactile frequency discrimination is not solely driven by mechanoreceptive afferents (Kuroki et al. 2017;Birznieks et al. 2019). There may be an additional system for vibrotactile frequency processing, possibly involving horizontal connections (Schwark and Jones 1989) or the secondary somatosensory cortex (Nelson et al. 2004;Chung et al. 2013;Kalberlah et al. 2013). Further research is needed to fully characterize S1 pRF properties as a function of frequency of vibrotactile stimulation.
In contrast to pRF size, we did not find that the amplitude of the BOLD signal was significantly affected by frequency of vibrotactile stimulation despite the substantial difference in kinetic energy delivered to cutaneous mechanoreceptors. Previous studies, however, reported that the BOLD amplitude can either increase (Nelson et al. 2004;Goloshevsky et al. 2008) or decrease (Chung et al. 2013) for increasing vibrotactile frequencies of stimulation. Especially when applying a vibrotactile stimulus for extended time periods, adaptation processes might have a negative effect on the BOLD signal amplitude. For the current experiments, we used an intermittent stimulation paradigm to minimize putative adaptation to the vibrotactile stimulus. It is possible that the current stimulation duration in combination with the intermittent stimulation paradigm equalized effects of different vibrotactile frequencies on BOLD amplitude.

Fingertip pRF size
We find that fingertip representations differ in pRF size. On average, cortical representations of the thumb exhibited the smallest pRF sizes, as we have reported previously (Schellekens et al. 2018). A gradual increase in pRF size is observed when progressing along the somatotopy, i.e. pRF size thumb < index < middle < ring < little finger. In a recent study, Puckett et al. (2020) reported larger pRF sizes in S1 for little finger representations compared to the index, middle and ring finger following a tactile stimulus, while measurements of the thumb were not included in their study. However, they did not observe a gradual change in pRF size across finger representations. The difference in results could possibly have been caused by methodological differences such as the smoothing applied in their analysis, which will generally increase pRF size estimates and increase the resemblance of pRF properties across voxels due to the Gaussian weighted average of neighboring voxels' timeseries in Gaussian smoothing algorithms. Additionally, the use of a separately estimated HRF in our study plausibly leads to better pRF estimations than using a canonical HRF as was done in the study of Puckett et al. (2020).
The difference in pRF size across fingertips occurred in all included Brodmann areas and under all vibrotactile frequency conditions. This makes it unlikely that the effect of fingertip representation on pRF size reflects functional hierarchical processes. Rather, the pRF size reflects the amount of integration of mechanoreceptive afferents from all fingers within single neuronal populations. Thus, the differences in pRF size per fingertip representation may be analogous to the increase in pRF size found in visual cortex for eccentricity representations, where foveal representations display smallest pRF sizes and outer eccentricities display larger pRF sizes (Smith et al. 2001;Dumoulin and Wandell 2008;Harvey and Dumoulin 2011). Assuming that neuronal populations representing the fovea might require high specificity for visual stimulus processing, a similar requirement may apply to somatosensory processing of tactile stimulation from the thumb and index finger. The thumb and index finger have the highest degree of motor acuity (Lachnit and Pieper 1990) and spatial acuity for somatosensory discrimination (Vega-Bermudez and Johnson 2001). Cortical pRF size might, additionally, relate to lower detection thresholds for thumb and index finger compared to the other digits in tactile discrimination tasks (Tamè et al. 2014). Our results indicate that neuronal populations that respond preferentially to the thumb and index finger receive relatively less mechanoreceptive input from the other fingers, compared to the cortical middle, ring and little finger representations, respectively.

Conclusion
We applied pRF modeling to investigate hierarchical information processing in S1 following vibrotactile stimulation of the five fingertips. PRF modeling allows for the assessment of a fingertip somatotopy in Brodmann areas BA3, BA1, and BA2. The pRF size portrays the degree of spatial information integration from the five fingertips within neuronal populations of cyto-architecturally distinct areas; smaller pRFs are associated with less spatial integration and earlier stages of the cortical processing hierarchy. pRF sizes were smallest in BA3, slightly increased for BA1, and approximately doubled in BA2, consistently across three different vibration frequencies. Additionally, we observed a difference in the time course of the hemodynamic response function among these Brodmann areas, with the shortest time-to-peak in BA3. Our findings confirm that the cortical hierarchy of the separate Brodmann areas in human S1 resembles the processing order observed in animal studies progressing from BA3 to BA1 and finally BA2, independent of the activated mechanoreceptors. Data availability All data can be made available.

Conflict of interest There are no conflicts of interest.
Ethical approval This study was approved by the local medical ethics committee.
Consent to participate All participants gave written informed consent prior to inclusion.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.