Fast frequency modulation is encoded according to the listener expectations in the human subcortical auditory pathway

The stimuli were three fast FM-sweeps: One sweep with a fast negative FM-rate (frequency span $\Delta f=-200$ Hz), one with a fast positive FM-rate (frequency span $\Delta f=200$ Hz), and one with a slow positive FM-rate (frequency span $\Delta f=100$ Hz; Fig 1A-B). We used 50 ms long sweeps in the frequency range of $f\sim 1500\,$Hz so that they had the typical properties of formant transitions in speech [47]. The sweep average frequencies were adjusted so that all FM-sweeps elicited the same average activity along the tonotopic axis and were perceived as having the same pitch [58, 57]; this design guaranteed that FM-direction and FM-rate selective neurons were necessary to differentiate between any two sweeps in the paradigm.

We arranged the stimuli in sequences of 8 FM-sweeps with 7 repetitions of the same sweep (standard) and one deviating sweep (deviant) (Fig 1C). Participants were instructed to report, with a button press, the position of the deviant within the sequence as fast and accurately as possible after identifying the deviant. Each sequence was characterised by the position of the deviant and $\Delta^{2}f=|\Delta f_{\text{deviant}}-\Delta f_{\text{standard}}|$: the absolute difference between the frequency spans of the deviant and the standard. With the three FM-sweeps we built three different combinations: one with $\Delta^{2}f=100$ Hz, where the two FM-sweeps have the same direction (up) but different modulation rate; one with $\Delta^{2}f=300$ Hz, where the two FM-sweeps have different rate and different direction; and one with $\Delta^{2}f=400\,$Hz, where to two FM-sweeps have the same modulation rate ($|\Delta f|=200$ Hz) but different direction (Fig 1B).

Expectations for each of the deviant positions were modulated by two abstract rules that were disclosed to the participants: 1) all sequences have a deviant, and 2) the deviant is always located in position 4, 5, or 6. The three deviant positions were used the same number of times along the experiment, so that the three deviant positions were equally likely at the beginning of the sequence. Therefore, the likelihood of finding a deviant in position 4 ($dev4$) after hearing 3 standards is $1/3$. However, if the deviant is not located in position 4, it must be located in either position 5 or 6, which makes the likelihood of finding a deviant in position 5 ($dev5$) after hearing 4 standards $1/2$. The likelihood of finding a deviant in position 6 ($dev6$) after hearing 5 standards is $1$.

To address the first research question, whether neural populations of human IC and MGB encode FM-rate and FM-direction, we tested whether these two nuclei show SSA to the FM-sweeps used in the experiment; namely, if neural responses in IC and MGB adapt to repeated FM-sweeps while preserving high responsiveness to FM-sweeps that deviate from the standards in FM-rate or FM-direction (Fig 1D). Since all sweeps were designed to elicit the same average activation across the tonotopic axis and elicited the same pitch percept, neural populations showing SSA to these FM-sweeps necessarily comprise neurons that are sensitive to FM-rate and FM-direction.

To address the second research question, whether IC and MGB responses encode FM-sweeps as prediction error with respect to the listener expectations, we used Bayesian model comparison. We considered two models. The first model assumed that adaptation to repeated fast FM-sweeps was mostly driven by habituation to the stimulus sequence properties, independently of participant’s expectations; namely, that neural populations habituate to repetitions of the standard, but show recovered responses to deviant irregardless of their position (habituation hypothesis; Fig 1E, h1). The second model assumed that adaptation was driven by predictive coding; namely, that neural responses to the deviants reflect prediction error with respect to the expectations of the participants (predictive coding hypothesis; Fig 1E, h2).

We measured BOLD responses in participants’ IC and MGB with an fMRI-sequence at 3-Tesla. The sequence was optimised to measure BOLD at relatively high spatial resolution (1.75 mm isotropic) while maintaining a high SNR (around 25).

All participants showed accuracies over 0.96 to all deviant positions (Fig 2A). Accuracy was slightly higher for the two more expected deviant positions, but differences between conditions were not significant ($p>0.1$, uncorrected). Reaction times (Fig 2B) showed a behavioural benefit of expectations: Participants reacted faster to more expected deviants (average $RT=770$ ms, $558$ ms and $246$ ms for deviants at positions 4, 5, and 6, respectively (Fig 2C); all differences were significant with $p<0.0001$, corrected for 3 comparisons).

We first studied whether the IC and MGB show SSA to fast FM-sweeps to test if the two nuclei are sensitive to FM-rate and FM-direction in humans. We estimated BOLD responses using a general linear model (GLM) with 6 different regressors: the first standard ($std0$), the standards after the first standard but before the deviant ($std1$), the standards after the deviant ($std2$), and deviants at positions 4, 5, and 6 ($dev4$, $dev5$, and $dev6$, respectively; Fig 1D). The conditions $std1$ and $std2$ were parametrically modulated [59] according to their positions to account for possible variations in the responses over subsequent repetitions (see Methods and Fig S1).

To compute SSA, we determined which voxels within the ICs and MGBs adapted to the standard (i.e., adaptation) and recovered responsiveness to deviants (i.e., deviant detection). SSA regions were then defined as the intersection between adaptation and deviant detection regions. ICs and MGBs were identified based on structural MRI data and an independent functional localiser (see Methods; IC and MGB ROIs; coloured patches in Fig 3). Within these ROIs, we used non-parametric ranksum tests ($N=18$; one sample per participant) to find which voxels showed significant adaptation to repeated standards (contrast $std0>0.5\,std1+0.5\,std2$). The associated $p$-maps were thresholded so that the false-discovery-rate $\text{\emph{FDR}}<0.05$. Surviving voxels constituted the adaptation ROIs (blue and purple patches in Fig 3). The same procedure was used to delimit the deviant detection ROIs (red and purple patches in Fig 3): the set of voxels within each anatomical ROI that responded significantly stronger to deviants than to repeated standards (contrast $dev4>0.5\,std1+0.5\,std2$; note that we compare the responses to the repeated standards with $dev4$ as this is the deviant position for which participants have the lowest expectation). The four anatomical ROIs showed significant adaptation (peak $p\leq 0.0001$) and deviant detection (peak $p<0.0001$; cluster size, exact peak $p$-values and MNI coordinates are shown in Table 1; all $p$-values corrected for four comparisons).

SSA regions were computed combining the unthresholded adaptation and deviant detection $p$-maps. The uncorrected $p$-value associated to SSA for a given voxel was $p_{\text{SSA}}=\max{(p_{\text{adaptation}},p_{\text{deviant detection}})}$. SSA $p$-maps where thresholded to $FDR<0.05$ to compute the SSA ROIs (Fig 3, purple). The four anatomical ROIs had extensive SSA regions (cluster sizes larger than 90 $mm^{3}$; peak $p\leq 0.0003$; exact peak $p$-values and MNI coordinates are shown in Table 1; all $p$-values corrected for four comparisons).

Significant SSA was also found at the single-subject level in 15 of the 18 participants ($p\leq 0.048$ for each of the 15 participants, corrected for the 596 voxels included in a global subcortical auditory ROI that comprised bilateral IC and MGB), but not all participants showed significant SSA in all ROIs (IC-L: 8 participants, $p\leq 0.049$; IC-R, MGB-L, MGB-R: 6 participants each, with $p\leq 0.048$; all $p$-values corrected for the number of voxels in the ROI and further corrected for four ROIs).

In the next step, we specifically tested whether the IC and MGB are similarly sensitive to FM-rate and FM-direction. To do that we analysed the data corresponding to: 1) trials where the standard and deviant differed only in modulation direction but not in absolute modulation rate; and 2) trials where the standard and deviant differed only in modulation rate but not in direction. If IC and MGB encode direction and rate, we would expect similar results in both partitions of the data. Conversely, if human IC and MGB are only sensitive to one of the two properties, we would expect null effects in the partition of the data where the standard and deviants differ only in that property.

Results were similar in both partitions of the data (Fig 4), demonstrating that the human IC and MGB encode both FM-direction and FM-rate.

We further corroborated that the levels of SSA were comparable for both types of FM changes at the single-subject level. In order to characterise FM-sensitivity with a number for each subject and FM-sweep combination, we used the SSA index [30] $SI$ (Eq (1); note that $SI>0$ is equivalent to the deviant detection contrast used in Fig 3).

We measured the difference in $SI$ to FM-direction ($SI_{dir}$) and FM-rate ($SI_{rate}$) in the voxels of the subject-specific SSA regions calculated in the previous section for each of the 15 subjects for which we obtained significant SSA. If FM-direction and FM-rate are both encoded in IC and MGB, we would expect no difference between these two partitions of the data. We measured the difference using Cohen’s $d=(\langle SI_{dir}\rangle-\langle SI_{rate}\rangle)/\sigma$), where $\langle SI\rangle$ is the average of $SI$ and $\sigma$ is the pooled standard deviation. The difference ranged between $d\geq-0.33$ and $d\leq 0.475$ across participants. The expected value of the difference ($E[d]=0.02\pm 0.05$) overlapped with zero, indicating once again that both FM-direction and FM-rate are already encoded in the subcortical auditory pathway.

To address our second question, we evaluated whether the average BOLD responses to deviants in the three different positions were affected by participant’s subjective expectations within the SSA regions. In congruence with the predictive coding hypothesis (Fig 1E, h2), the response profile showed reduced responses for more expected deviants (Fig 5).

Formal statistical testing confirmed that responses to different deviant positions were different in all ROIs for all contrasts among deviant positions: $dev4\neq dev5$ ($|d|\geq 0.99$ and $p<0.006$), $dev4\neq dev6$ ($|d|\geq 2.39$ and $p<0.00005$), and $dev5\neq dev6$ ($|d|\geq 1.74$ and $p<0.0003$; all $p$-values corrected for $3\times 4$ 12 comparisons). Exact $p$-values and effect sizes are listed in Table 2. All statistical tests included one sample per participant, ROI, and deviant position.

To corroborate that differences were present at the single-subject level we run a correlation analysis for each of the 15 participants for which we obtained significant SSA. In each participant, we computed the Pearson’s correlation between the BOLD responses elicited by each deviant location with its likelihood of occurrence (namely, $1/3$ for deviant 4, $1/2$ for deviant 5, and $1$ for deviant 6). If BOLD responses reflect prediction error, we would expect a negative correlation between the likelihood and the responses. We found significantly negative correlations in all 15 participants ($\rho\in[-0.87,-0.42]$, all $p<0.03$; all Pearson tests had $9\times 3=27$ samples, 3 per run).

We used Bayesian model comparison to formally evaluate the response properties in each voxel of the IC and MGB ROIs. This approach provides for a quantitative assessment of the likelihood that each of the two hypotheses (Fig 1E) can explain the responses in each voxel. This analysis is sensitive to possible region-specific effects that could have been averaged out when aggregating the z-scores across voxels in each ROI.

Following the methodology described in [60, 61], we first calculated the log-likelihood of each model in each voxel of the two ICs and MGBs in each participant. Each model yields different predictions on the relative amplitudes to different positions in the sequences (Fig 1E). We tested h1 and h2 to adjudicate between the habituation and predictive coding explanations of the responses. H1 assumed an asymptotic decay of the standards and recovered responses to the deviants; h2 assumed that the responses to both deviants and standards would depend on the participant’s expectations (Fig 6; for exact values, see Methods). Participant-specific log-likelihoods were used to compute the Bayes factor $K$ (i.e., the ratio of the posterior likelihoods) between h1 and h2.

H2 was the best explanation for the data in the majority of voxels of the four ROIs (Figures 7 and 8): h2 was more likely than h1 in 99% and 80% of the left and right IC, respectively, and in all voxels of the left and right MGB.

To test whether the effect was present at the single-subject level, we computed $K$ independently for each subject in the subject-specific bilateral IC and MGB. Since we performed the group analyses over the entire ROIs (and not only the SSA regions), here we used the full anatomical ROIs of each participant. We measured for how many voxels within each region and participant h2 was the better explanation of the data. H2 was the better explanation of the data (more than 50% of voxels) in 16 (IC-L), 13 (IC-R), 13 (MGB-L), and 15 (MGB-R) participants. In 7 (IC-L), 5 (IC-R), 10 (MGB-L), and 8 (MGB-R) participants H2 was the better explanation even in more than 75% of voxels.

The auditory pathway is anatomically subdivided into two sections: the primary (lemniscal) or secondary (non-lemniscal) pathways. The primary pathway is characterised by neurons that carry auditory information with high fidelity and it is generally regarded as responsible for the transmission of bottom-up sensory input [48]. The secondary pathway has wider tuning curves and it is generally regarded as responsible for the integration of contextual and multisensory information [48].

Both IC and MGB comprise regions that participate in both, the primary and secondary pathways [48]. The primary subdivision of the IC is its central nucleus, while the cortices constitute the secondary subdivisions. The primary subdivision of the MGB is its ventral section, while the medial and dorsal sections constitute the secondary subdivisions.

In rodents, SSA and prediction error to pure tones are significantly stronger in secondary subdivisions (e.g., [31]). In humans, prediction error is similarly strong in primary and secondary MGB for pure tones [40]. Here we test for differential representations of prediction error to FM-sweeps in MGB.

Distinguishing between the primary and secondary subsection of the IC and MGB non-invasively is technically challenging [62]. A recent study [63] distinguished two distinct tonotopic gradients of the MGB. The ventral tonotopic gradient was identified as the ventral or primary (vMGB) subsection of the MGB (see Fig 9A, green). Although the parcellation is based only on the topography of the tonotopic axes and their anatomical location, the region is the best approximation to-date of the vMGB in humans. No parcellation of the IC is available to-date.

Both primary and secondary subdivisions of bilateral MGB showed SSA. SSA strength was measured in each voxel using the SSA index (Eq (1)). Distributions of the $SI$ across the voxels of each of the subdivisions were comparable in both hemispheres (Fig 9B) demonstrating that SSA is not confined to nor stronger in the secondary MGB.

Predictive coding (h2) was the best explanation for the responses of all voxels in the two subdivisions of the left MGB, and in 95% and 97% of the primary and secondary subdivisions of the right MGB, demonstrating that FM-sweeps are encoded as prediction error in both, primary and secondary subdivisions of bilateral MGB. Moreover, the distributions of the Bayes’ factor $K$ between the predictive coding (h2) and adaptation (h1) hypotheses were comparable across subdivisions (Fig 9C).

To study whether the same neural populations are in charge of encoding prediction error to FM-sweeps and pure tones, we compared the topographic distribution of the Bayes factor $K$ between the h2 and h1 in our data with the topographic distribution of the Bayes factor $K$ we obtained in a previous experiment, where we measured BOLD responses to the same experimental paradigm as here but using pure tones [40]. We computed the correlation between both $K$ across voxels of each of the four ROIs, as defined by the anatomical atlas from [64].

Distribution of $K$ to both families of stimuli was strongly correlated across voxels of bilateral IC (left, $\rho=0.47$, $p=10^{-7}$; right, $\rho=0.34$, $p=10^{-4}$; $p$-values corrected for 4 comparisons), but not across voxels of the MGBs (left, $\rho=-0.11$, $p=0.22$; right, $\rho=0.15$, $p=0.1$; uncorrected $p$-values).

Eighteen German native speakers (12 female), aged 19 to 31 years (mean 24.6), participated in the study. None of them reported a history of psychiatric or neurological disorders, hearing difficulties, or current use of psychoactive medications. Normal hearing abilities were confirmed with pure tone audiometry (250 Hz to 8000 Hz); all participants had hearing threshold equal to or below 15 dB SPL in the frequency range of the FM sweeps (1000 Hz-3000 Hz). Participants were also screened for dyslexia (German SLRT-II test [97], RST-ARR [98], and rapid automatised naming (RAN) test of letters, numbers, objects, and colours [99]) and autism (Autism Spectrum Quotient; AQ [100]). All scores were within the neurotypical range (SLRT: $\min(\max(PR_{\text{words}},PR_{\text{pseudowords}}))=21$, higher than the cut-off value of 16, following the same guidelines as [101]; RST-ARR: all $PR\geq 31$, higher than the cut-off value of 16; RAN: maximum of 3 errors and $RT<36$ seconds in each of the four categories; AQ: all participants $AQ\leq 31$, under or equal to the cut-off value of 32).

Since we had no estimations of the possible sizes of the effects, we maximised our statistical power by recruiting as many participant as we could fit in the MRI measurement time allocated to the study. This number was fixed to twenty before we started data collection, but two participants dropped out of the study during data collection. We maximised the amount of data collected for participant to reduce random error to a minimum and maximise the likelihood of measuring effects at the single-subject level.

All sounds were 50 ms long (including 5 ms in/out ramps) sinusoidal FM sweeps. The frequency sweeps lasted for 40 ms and were preceded and followed by 5 ms long segments of constant frequency that overlapped with the in/our ramps. The sweeps were assembled in the frequency space to avoid discontinuities of the final waveforms.

We used a total of three sweeps during the experiment: a fast up sweep with starting frequency $f_{0}=1000$ Hz and ending frequency $f_{1}=1200\,$Hz ($\Delta f=200$ Hz); a slow up sweep with $f_{0}=1070$ Hz and $f_{1}=1170$ Hz ($\Delta f=100$ Hz), and a fast down sweep with $f_{0}=1280$ Hz and $f_{1}=1080$ Hz ($\Delta f=-200$ Hz). The sweeps had different average frequencies to ensure that they elicited the same average spectral activity along the tonotopic axis and the same pitch percept (see [58] for details).

From those 3 sweeps we constructed 6 standard-deviant frequency combinations that were used the same number of times across each run, so that all sweeps were used the same number of times as deviant and standards. Each sequence consisted of 7 repetitions of the standard sweep and a single instance of the deviant sweep. Sweeps were separated by 700 ms inter-stimulus-intervals (ISI), the shortest possible ISI that allowed the participants to predict the fully expected deviant [40], amounting to a total duration of 5300 ms per sequence.

In each trial of the fMRI experiment, participants listened to one tone sequence and reported, as fast and accurately as possible using a button box with three buttons, the position of the deviant (4, 5 or 6). The inter-trial-interval (ITI) was jittered so that deviants were separated by an average of 5 seconds, up to a maximum of 11 seconds, with a minimum ITI of 1500 ms. We chose such ITI properties to maximise the efficiency of the response estimation of the deviants [102] while keeping a sufficiently long ITI to ensure that the sequences belonging to separate trials were not confounded.

All but one participant completed 9 runs of the main experiment across three sessions; participant 18 completed only 8 runs for technical reasons. Each run contained 6 blocks of 10 trials. The 10 trials in each block used one of the 6 possible sweep combinations, so that all the sequences within each block had the same standard and deviant. Thus, within a block only the position of the deviant was unknown, while the deviant’s FM-direction and FM-rate were known. The order of the blocks within the experiment was randomised. The position of the deviant was pseudorandomised across all trials in each run so that each deviant position happened exactly 20 times per run but an unknown amount of times per block. This constraint allowed us to keep the same a priori probability for all deviant positions in each block. In addition, there were 23 silent gaps of 5300 ms duration (i.e., null events of the same duration as the tone sequences) randomly located in each run [102]. Each run lasted around 10 minutes, depending on the reaction times of the participant.

Due to an undetected bug in the presentation code, information on the exact sweep combination used in each trial was unavailable for some runs. The bug affected the first three runs of participants 1, 2, 4, and 5; and the first six runs of participant 3. This information was not relevant for the analyses that aggregated the data across sweep combinations, and affected only the analyses of Fig 4, where we excluded the affected runs of participants 1, 2, 4, and 5, and participant 3 altogether.

We also run a functional localiser that was designed to activate the participant’s IC and MGB. Each run of the functional localiser consisted on 20 blocks of 16 seconds and lasted for about 6.5 minutes. Ten of the blocks were silent; the remaining blocks consisted on presentations of 16 sounds of one second duration each. Sounds were taken from a collection of 85 natural sounds collected by [62]. Participants were instructed to press a key when the same sound was repeated twice to ensure that they attended the sounds; behavioural data from the functional localiser was not used in the analysis.

Each session consisted on three runs of the main experiment, interspersed with two runs of the functional localiser. All runs were separated by breaks of a minimum of 1 minute to allow the participants rest. Fieldmaps and a whole-head EPI (see Data acquisition) were acquired between the third and fourth run. In the first session, we also measured an structural image before the fieldmaps. The first run of the first session was preceded by a practice run of four randomly chosen trials to ensure the participants had understood the task. We acquired fMRI during the practice run in order to allow the participants to undertake the training with MRI-noise.

MRI data were acquired using a Siemens Trio 3 T scanner (Siemens Healthineers, Erlangen, Germany) with a 32-channel head coil.

Functional MRI data were acquired using echo planar imaging (EPI) sequences. We used partial coverage with 24 slices. The volume was oriented in parallel to the superior temporal gyrus such that the slices encompassed the IC, the MGB, and the superior temporal gyrus. In addition, we acquired one volume of an additional whole-head EPI with the same parameters (including the FoV) and 84 slices during resting to aid the coregistration process (see Data preprocessing).

The EPI sequence had the following acquisition parameters: TR = 1900 ms, TE = 42 ms, flip angle 66^∘, matrix size $88\times 88$, FoV 154 mm$\times$154 mm, voxel size 1.75 mm isotropic, bandwidth per pixel $1386$ Hz/px, and interleaved acquisition. During functional MRI data acquisition, cardiac signal was acquired using a scanner pulse oximeter (Siemens Healthineers, Erlangen, Germany).

Structural images were recorded using an MPRAGE [103] T1 protocol with 1 mm isotropic resolution, TE = 1.95 ms, TR = 1000 ms, TI = 880 ms, flip angle 1 = 8^∘, FoV = 256 mm$\times$256 mm.

Stimuli were presented using MATLAB (The Mathworks Inc., Natick, MA, USA) with the Psychophysics Toolbox extensions [104] and delivered through an Optoacoustics (Optoacoustics Ltd, Or Yehuda, Israel) amplifier and headphones equipped with active noise-cancellation. Loudness was adjusted independently for each participant to a comfortable level before starting the data acquisition.

The preprocessing pipeline was coded in Nipype 1.5.0 [105], and carried out using tools of the Statistical Parametric Mapping toolbox, version 12; Freesurfer, version 6 [106]; the FMRIB Software Library, version 5 (FSL) [107]); and the Advanced Normalization Tools, version 2.3 (ANTS) [108]. All data were coregistered to the Montreal Neurological Institute (MNI) MNI152 1 mm isotropic symmetric template.

First, we realigned the functional runs. We used SPM’s FieldMap Toolbox to calculate the geometric distortions caused in the EPI images due to field inhomogeneities. Next, we used SPM’s Realign and Unwarp to perform motion and distortion correction on the functional data. Motion artefacts, recorded using SPM’s ArtifactDetect, were later added to the design matrix (see Estimation of the BOLD responses).

Next, we used Freesurfer’s recon-all routine to calculate the boundaries between grey and white matter (these are necessary to register the functional data to the structural images) and ANTs to compute the transformation between the structural images and the MNI152 symmetric template.

Last, we coregistered the functional data to the structural image with Freesurfer’s BBregister, using the boundaries between grey and white matter of the structural data and the whole-brain EPI as an intermediate step. Data was analysed in the participant space, and then coregistred to the MNI152 template. Note that, since the resolution of the MNI space (1 mm isotropic) was higher than the resolution of the functional data (1.75 mm isotropic), the transformation resulted in a spatial oversampling.

All the preprocessing parameters, including the smoothing kernel size, were fixed before we started fitting the general linear model (GLM) and remained unchanged during the subsequent steps of the data analysis.

Physiological (heart rate) data was processed by the PhysIO Toolbox [109], that computes the Fourier expansion of each component along time and adds the coefficients as covariates of no interests in the model’s design matrix.

First level analyses were coded in Nipype and carried out using SPM. Second level analyses were carried out using custom code in MATLAB. The coregistered data were first smoothed using a 2 mm full-width half-maximum kernel Gaussian kernel with SPM’s Smooth.

The first level GLM’s design matrix for the main experiment included 6 regressors: first standard (std0), standards before the deviant (std1), standards after the deviant (std2), and deviants in positions 4, 5, and 6 (dev4, dev5, and dev6, respectively; Fig 1). Conditions std1 and std2 were modelled using linear parametric modulation [59], whose linear factors were coded according to the position of the sound within the sequence to account for effects of habituation [40]. The first level GLM’s design matrix for the functional localiser included 2 conditions: sound and silence. On top of the main regressors, the design matrix also included the physiological PhysIO and artefact regressors of no-interest.

We used a recent anatomical atlas of the subcortical auditory pathway [64] to compute prior regions corresponding to the left IC, right IC, left MGB, and right MGB, respectively. The atlas comprises three different definitions of the ROIs calculated using 1) data from the big brain project, 2) postmortem data, and 3) fMRI in vivo-data. To compute the prior coarse region for each nuclei we combined the three masks and inflated the resulting regions with a Gaussian kernel of 1 mm fwdh.

The final IC and MGB regions were computed by combining the prior coarse regions with the results from the contrast $sound>silence$ of the functional localiser. Within each region, we thresholded the contrast to increasingly higher values until the number of surviving voxels equalled the volume of the region reported in [64]; namely, 146 voxels for each of the ICs, and 152 for each of the MGBs.

We used the coefficients of the GLM or beta estimates from the first level analysis to calculate the adaptation (Fig 3, blue patches) and deviant detection (red patches) ROIs, defined as the sets of voxels within the IC and MGB ROIs that responded significantly to the contrasts $std0>0.5std1+0.5std2$ and $dev4>0.5std1+0.5std2$, respectively. Significance was defined as $p<0.05$, family-wise-error (FDR)-corrected for the number of voxels within each of the IC/MGB ROIs. SSA voxels are defined as voxels that show both, adaptation and deviant detection; thus, we calculated an upper bound of the $p$-value maps for the SSA contrast as the maximum of the uncorrected $p$-values associated to the adaptation and deviant detection contrasts. The SSA ROIs (Fig 3, purple patches) were calculated by FDR-correcting and thresholding the resulting $p$-maps at $\alpha=0.05$. All calculations were performed using custom-made scripts (see Data and code availability).

The Bayesian analysis of the data consisted as well of first and second level analyses. In the first level, we used SPM via nipype to compute the log-evidence in each voxel of each participant for each of the four models (see Fig 6). The models were described using one regressor with parametric modulation whose coefficients corresponded to a simplified view of the expected responses according to each model (Table 3). The expected responses of each model were the same in all trials that had the same standard-deviant combination and deviant position. Given the model amplitude(s) $a_{n}$ and the timecourse of a voxel $y$, SPM calculates the log-evidence of the linear model $y=\sum\beta_{n}a_{n}+\xi$, where $\beta_{n}$ are the linear coefficients of each regressor and $\xi$ are noise terms.

Log-evidence maps for each participant were combined following the random-effects-equivalent procedure described in [60, 61] to compute the posterior probability maps associated to each model at the group level. We combined the maps using custom scripts (see Data and code availability). Histograms shown in Figures 7 and 9 are kernel-density estimates computed with the distribution of the posterior probabilities across voxels for each of the SSA ROIs.

All pairwise comparisons reported in the study were evaluated for significance using two-tailed Ranksum tests. Unless stated otherwise, $p$-values for all analyses that comprised multiple testing were corrected using the Holm-Bonferroni method. A result was deemed statistically significant when the corrected $p<0.05$.

Derivatives (beta maps and log-likelihood maps, computed with SPM) and all code used for data processing and analysis are publicly available in https://osf.io/f5tsy/.