{'fraction_non_alphanumeric': 0.03519369009584664, 'fraction_numerical': 0.015924520766773163, 'mean_word_length': 4.781529581529582, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spatial organization of frequency preference and selectivity in the human inferior colliculus

To date, the functional organization of human auditory sub-cortical structures can only be inferred from animal models. Here we use high-resolution functional MRI at ultra-high magnetic fields (7 Tesla) to map the organization of spectral responses in the human inferior colliculus (hIC), a sub-cortical structure fundamental for sound processing. We reveal a tonotopic map with a spatial gradient of preferred frequencies approximately oriented from dorso-lateral (low frequencies) to ventro-medial (high frequencies) locations. Furthermore, we observe a spatial organization of spectral selectivity (tuning) of fMRI responses in the hIC. Along isofrequency contours, fMRI-tuning is narrowest in central locations and broadest in the surrounding regions. Finally, by comparing sub-cortical and cortical auditory areas we show that fMRI-tuning is narrower in hIC than on the cortical surface. Our findings pave the way to non-invasive investigations of sound processing in human sub-cortical nuclei and to studying the interplay between sub-cortical and cortical neuronal populations.


Supplementary Figure S3. Group level cortical and sub cortical activation maps
The left column shows the group functional activation (F-map) in response to simple tones (first row), natural sounds (second row) and all stimuli together (third row) as obtained by the cortical analysis (1 mm isotropic voxels; FDR q = 0.05; data from Experiment 1) sampled on the average (after cortex based alignment) surface of the five subjects. The right column shows the group functional activation (F-map) in response to simple tones (first row), natural sounds (second row) and all stimuli together (third row) as obtained by the sub-cortical analysis (0.5 mm isotropic voxels; FDR q = 0.001) overlaid to the average (after Talairach alignment) anatomy of the five subjects. Parcellation of auditory cortical areas based on group tonotopy and tuning width maps. Group functional maps from the natural sounds experiment of spectral preference (i.e. tonotopy; red-yellow-green-blue scale, the same map is repeated in the second and third row) and spectral tuning (i.e. bandwidth; purple-blue-greenyellow scale, fourth row) are sampled on the average (after cortex based alignment) cortical surface of the five subjects. In the second row, the white circle outlines the high-low-high main frequency gradients consistently observed across experiments. In the third and fourth row, full and dashed lines outline the distinction between primary auditory cortex and lateral belt areas. Numbering within the primary auditory cortex highlights the separation between hA1 and hR. Results of the single subject permutation analysis (N = 500) for the assessment of significance of center frequency (left column) and tuning width (right column) in the natural sounds modeling. Individual maps are thresholded at p < 0.01 uncorrected.

Supplementary Figure S9. Single subject cross validation analysis
Results of the single subject cross-validation analysis (20 splits) for the assessment of consistency in the assignment of best frequency (data of Experiment 2). Individual maps are thresholded at 60% (12 out of 20 splits) consistency.
Each condition was presented twice per run, and we acquired six runs per subject.
In Experiment 2, stimuli were derived from recordings of natural sounds (n = 144; e.g. speech sounds, vocalizations, animal cries, tools, nature sounds), and cut to a one second duration. Because the estimation of the voxels response profiles (see below) requires measuring the responses to a large number of stimuli, in this experiment stimuli were presented according to a fast event-related scheme. Despite the possible influence of acoustic scanner noise, we selected such a scheme -and not e.g. a "sparse" sampling approach 49 -because of its high temporal efficiency (which is crucial when the presentation of many stimuli is required) and based on our previous experience with cortical tonotopic maps at 3 Tesla 31 . Note that to ensure a clear perception of the sounds during scanning, in this experiment (as well as in Experiment 1 and 3) all stimuli were presented in silent gaps in between the acquisitions of two subsequent brain volumes 38, 50,51 Each stimulus was presented three times across six different runs, and the inter-stimulus-interval (ISI) was jittered between 5200 and 10400 ms. Zero trials (where no sound was presented, 6% of total trials) and catch trials (repeating the sound of the previous trial, 5% of total trials) were added. Subjects were asked to press a button when a sound was repeated; these catch trials were excluded from the analysis of the data. Five subjects (median age = 32, 3 males) participated in both Experiment 1 and 2.
In Experiment 3, 800 ms long amplitude-modulated (8 Hz; Each block consisted of six sounds, lasted 18 seconds and was separated by 12 seconds of no stimulation from the following block. Each condition was presented twice per run, and we acquired six runs per subject. Four subjects (median age = 30; 3 males) were scanned in Experiment 3.
In all three experiments, sounds were sampled at 44.1 KHz. Sound onset and offset were ramped with a 10 ms linear slope and sound energy (RMS) was equalized.
Prior to scanning, all sounds were subjectively equalized for loudness. In Experiments 1 and 2 sounds were presented with an MR compatible audio system based on air tubes (Avotec Inc.) with a linear frequency transfer function up to about 4 KHz. The system was modified in such a way that the final piece of the air tubes was directly inserted n the subjects' ear through moldable earplugs. In experiment 3 sounds were presented with an MR compatible piezoelectric/earbud system (Sensimetric Inc.) with a linear frequency transfer function up to about 8 KHz. In all experiments sounds were presented at about 60 dB. All subjects were asked to place memory foam pieces cut to an average ear size on top of their ears for additional sound protection.

Imaging parameters
Data were acquired on a 7T whole body system driven by a Siemens console using a head gradient insert operating at up to 80 mT/m with a slew rate of 333 T/m/s. A head RF coil (single transmit, 16 receive channels) was used to acquire anatomical (T 1 weighted) and functional (T 2 * weighted BOLD) images. In each experiment, T 1 weighted (1 mm 3 ) images were acquired using a modified Magnetization Prepared Rapid Acquisition Gradient Echo (MPRAGE) sequence (TR=2500 ms; TI=1500 ms; TE=3.67 ms). Proton density (PD) images were also acquired together with the T 1 weighted images (both acquisitions are interleaved in the modified MPRAGE sequence) and were used to minimize inhomogeneities in T 1 weighted images 52 .
Acquisition time for anatomical images was ~ 7 minutes. T 2 * weighted functional data were acquired using a clustered Echo Planar Imaging (EPI) technique. In Experiment 1 and 3, each brain volume was scanned at 1.5 mm isotropic resolution (44 slices; GRAPPA acceleration X3; partial Fourier 6/8; TR = 3000 ms; time of acquisition [TA] = 1500 ms). In Experiment 2, brain volumes were scanned at 1.5 mm isotropic resolution (31 slices; GRAPPA acceleration X3; partial Fourier 6/8; TR = 2600 ms; time of acquisition [TA] = 1200 ms). For all experiments, slices were placed in order to maximize efficiency (short acquisition time) and cover the brain transversally from the inferior portion of the anterior temporal pole to the superior portion of the superior temporal gyrus bilaterally (see Supplementary Fig. 1). All functional sessions consisted of six runs of approximately 10 minutes each. As cardiac and respiration cycles were not recorded, we did not apply any post-hoc correction to reduce cardiac induced pulsation noise that may decrease the temporal SNR in brainstem regions. Despite this we found reliable activation in the IC to all presented sounds (see below).
Functional and anatomical data were analyzed with BrainVoyager QX.
Preprocessing consisted of slice scan-time correction (with sinc interpolation), temporal high pass filtering (removing drifts of 3 cycles or less per run in data collected from Experiments 1 and 3 and of 7 cycles or less per run in data collected from Experiment 2) and 3-dimentional motion correction. Additional spatial smoothing (Gaussian kernel with FWHM = 2 mm) and temporal smoothing (2 consecutive data points) were applied to the data collected in Experiments 1 and 3.
Functional data were co-registered to the anatomical data collected in the same experiment. Anatomical and functional data were co-registered across sessions and normalized in Talairach space. Functional data were resampled (with sinc interpolation) in the normalized space at a resolution of 0.5 mm isotropic and, for the hIC analysis, limited to a bounding box (60 x 70 x 60 mm) surrounding the inferior colliculus bilaterally. Supplementary Fig. 6 and Supplementary Fig. 7 show the effects of the pre-processing stages (i.e. spatial smoothing and resampling) and map level spatial smoothing on the tonotopic maps obtained for a single subject in Experiment 3. While improving the visualization of the maps, the pre-processing strategies did not change relevantly the overall tonotopic spatial gradient.

Functional data analysis
The response to simple tones and natural sounds was evaluated by computing a fixed

Simple tones experiments.
The spatial layout of frequency selective responses elicited by simple tones (i.e. Experiments 1 and 3) was obtained by computing "best frequency" maps 13 . A single subject GLM analysis was used to compute the response to the three center frequencies (.5; .1.5; 2.5 kHz) in Experiment 1 and to the eight center frequencies (.18; .30; .51; .86; 1.46; 2.48; 4.19; 7.09 kHz) in Experiment 3. The voxels' best frequency (BF) was assigned for each run separately as the frequency to which the voxel responded strongest. A leave-one-run out cross-validation procedure was used to assess the reproducibility of these maps within each single subject. Individual subject tonotopy maps were created by averaging the voxels' BF across runs, and by color-coding the resulting map according to a linear red-yellow-green-blue color scale. Group maps were obtained by averaging single subject BF maps.

Natural sounds experiment.
Tonotopic and tuning width maps from Experiment 2 were obtained using an "encoding" technique 20,21,23 that allows estimating the voxels' response to features (e.g. frequency) from the fMRI patterns elicited by complex sounds. First, the frequency content of all natural sounds (n = 144) was characterized by means of a Fourier analysis. This resulted in a matrix W (n x f) of Fourier coefficients were f is the number of resulting frequency bins. Frequencies ranged linearly from 0 to 8 kHz and the number of frequency bins was fixed to 40. Second, the response pattern Y (v x n) to each natural sound across the v voxels in the fMRI data set was estimated using a GLM analysis. For each voxel j a response vector Y j (1 x n) was computed in two steps. First, the hemodynamic response function (HRF) common to all stimuli was obtained using deconvolution. Then, the voxel specific HRF was used to estimate the response to all individual sounds using a GLM analysis with one predictor per sound 53 . Finally, the matrix R (v x f) of voxels' frequency response profiles was obtained by estimating, for each voxel independently, the linear model that links the frequency content of the sounds W and the single voxel fMRI response to all sounds Y j : (1) where R j (1 x f) represents the contribution of each frequency bin to the overall response of voxel j (i.e. the voxels' spectral response profile [SRP]). The solution to this linear model was obtained using Relevance Vector Machine (RVM) regression 54 .
A 20-fold cross-validation procedure was used leaving out 24 sounds per fold.
Tonotopic and tuning width maps were obtained as the mean across the 20 folds.
Note that this analysis requires solving (per each voxel and fold) a linear system (Eq. As the tonotopy maps obtained from tones (Experiments 1 and 3) were estimated using a simple least squares solution (i.e. General Linear Model), statistical significance of the estimates could be obtained using standard parametric testing.
Conversely, the tonotopy maps obtained from natural sounds (i.e. the voxels' SRP; Experiment 2) were obtained using a regularized regression approach (RVM) to robustly deal with the partial collinearity of the design matrix (W). Regularized models, such as RVM, do not have a null-distribution and thus do not allow for a simple parametric test. Here, we have estimated the empirical null-distribution by randomly permuting (N = 500 per subject) the voxels' responses to all sounds (Y) and repeating the fitting procedure. Note that as the permutation was performed after the deconvolution analysis and prior to the fitting procedure, the serial autocorrelation present in the original fMRI time series does not represent an issue here. For each fold, best frequency (BF) and tuning width (TW) were estimated per voxel using a two-step procedure. First, we computed the average voxel response within three frequency windows (center: .5; .1.5; 2.5 KHz; width 600 Hz). The center frequency (in Hz) of the maximum responding window was assigned as the voxels' BF.
Tonotopic maps were obtained by color-coding the estimated BF using a linear redyellow-green-blue color scale. Values obtained from the permutation testing were used to evaluate the significance of the BF assignment and the results are reported in Supplementary Fig. 8 (first column). In Supplementary Fig. 9 we report the results of a cross-validation procedure (N=20) that highlights the stability of our frequency estimates. Second, we fitted a Gaussian function with mean constrained within the frequency window eliciting the most consistent maximum response across all folds.
The goodness of fit of the Gaussian was compared to the goodness of fit obtained by the same procedure on the 500 permuted results to obtain significant maps for the tuning width ( Supplementary Fig. 8, second column). For each voxel j, the full width at half maximum (FWHM) of the fitted Gaussian in Hz was used to compute the tuning width (TW) as: (2) We observed, in correspondence with previous studies in a variety of animals, a positive correlation between maps of TW and BF 27, 55 . The dependence of TW on BF was estimated with a smoothing spline fit. The residuals of this fit (TW res ) were used to create final maps of tuning width. Tuning width maps were obtained by colorcoding the estimated TW res using a linear yellow-green-blue-purple color scale. Group maps were obtained by averaging single subject BF and TW res maps. Tonotopic and tuning width maps were obtained following the same color-coding used for single subject maps (tonotopy: red-yellow-green-blue; tuning width: yellow-green-bluepurple).
The correspondence between tonotopic maps obtained from simple tones (Experiment 1) and natural sounds (Experiment 2) was evaluated by computing the correlation between these two tonotopic maps in each subject and between group maps. Significance was obtained with permutation testing (n=1000) while accounting for the spatial correlation present in the single maps (i.e. by applying the same smoothing procedure to the unpermuted and permuted maps). In four of the five subjects and for the group maps, we observed a significant correlation between the tonotopic maps obtained from the two experiments (p<0.01 for individual subjects; over 3602 voxels corresponding to a p<0.05 across group maps).
In order to examine the pattern of tuning width in the hIC for individual subjects and the group map, we performed an additional analysis on each of these maps (for Experiments 1 and 3). First, for each map we selected eight equally spaced sagittal slices (four slices in left and right hemisphere each; slices selected in center of hIC to minimize the variation in size of the maps). Second, we defined the isofrequency contour at 1.5 kHz (i.e. middle frequency) in each slice. Third, all isocountours were linearly interpolated to hundred dimensional vectors and aligned to each other. The average length of these iso-contours was approximately 8 mm.
Finally, we sampled the tuning width map along each contour (see Supplementary Finally, to examine whether the spatial distribution of tuning width in hIC could result from vasculature differences throughout hIC, we evaluated the correlation between tuning width and the overall response significance to natural sounds (F-Map; data from Experiment 2) or to the beta-map representing the response amplitude to natural sounds. We found no correlation between these measures, which ensures that the pattern observed in tuning width maps cannot be explained by hemodynamic effects.

Gradient analysis on best frequency maps.
In order to quantify the direction of the main frequency gradient, a gradient analysis was performed on the single subject and group tonotopic maps of the hIC obtained from all experiments. For each voxel, we computed -in subsequent sagittal slicesthe direction of increasing best frequency with respect to a transversal (i.e. horizontal) vector pointing in the anterior to posterior direction. To avoid nuisance effect from voxels lying at the border of the map, the direction was averaged across voxels within a region of interest obtained by eroding the original hIC mask (i.e. eliminating the outer most shell of the hIC).

Parcellation of cortical auditory fields
Group tonotopic and tuning width maps were used to delineate three cortical auditory fields (hA1, hR, lateral belt; see Supplementary Fig. 5). We defined the narrow region along HG as the auditory core. Based on observations of coinciding cytoarchitectonic borders and functional characteristics 29 , this region was divided into a medial and lateral region at the low frequency gradient reversal. The resulting regions each contained a complete frequency gradient and may reflect the human homologues of monkey A1 and R 28,29 . Laterally and immediately adjacent to the auditory core, occupying Heschl's sulcus (HS) and the anterior part of planum temporale, we identified an auditory belt area. This area is distinguished from the auditory core by its tuning width, which is consistently broader 30,40-43 .
The cortical regions of interest as drawn on the average cortical surface of all five subjects were projected back to the individual volumes. The distribution of voxels across tuning width values in the hIC, hA1, hR and lateral belt regions in left and right hemispheres is reported in Fig. 5 for each individual subject. All distributions were normalized for the total number of voxels in each area. Note that, as TW res values were computed by removing the relation between TW and BF separately for cortical and subcortical regions, we could not use these values for comparisons across regions. Instead, we used uncorrected TW for this part of the analysis.
While in most subjects the tuning width in cortical regions was limited to 0.6 octaves (ie. TW=2.5) or broader, in the IC several voxels (consistently across subjects and hemispheres) were more narrowly tuned (up to 0.45 octaves; i.e. TW=3.23). We tested, in four intervals equally spaced between the minimum and maximum observed tuning widths, for differences in number of voxels across the four regions (see Fig. 6).
We observed significant differences, with the lateral belt presenting more broadly tuned voxels and the IC more narrowly tuned voxels. We also tested (non parametric test Wilcoxon rank-sum test) for differences in mean tuning between areas and across subjects and observed, in both left and right hemisphere, narrower tuning in the IC compared to lateral belt areas (p = 0.05; n = 5) and between right IC and right hA1 (p = 0.05 ; n = 5).