Introduction

Anatomical connections of the human cortex differ across cortical depths to support feedforward or feedback information processing between cortical and subcortical regions [1]. Neurophysiological evidence from invasive recordings suggests that the feedforward and feedback processing are differentially manifested in neuronal oscillations [2–5] and broadband high frequency (>80 Hz) neural activities (HFA) [6,7]. The link between structural and neural oscillatory features in feedforward and feedforward processing is that dominant neural oscillatory frequency changes between alpha, beta, and gamma bands across cortical depths [8–12]. Distinct early and late components of HFA at deep, intermediate, and superficial cortical layers were recently reported to support the feedforward and feedback auditory and visual processing in non-human primate [7]. However, the link between cytoarchitecture and neural oscillations as well as HFA to dissociate feedforward and feedback modulations has been less studied in humans.

This challenge can be mitigated by magnetic resonance imaging (MRI) for its non-invasive nature and submillimeter resolution. In particular, the intrinsic hemodynamic contrast in functional MRI (MRI) is closely related to local field potentials [13–15]. Previous fMRI studies in the visual and sensorimotor cortices reveal that fMRI signals are negatively correlated with alpha and beta band oscillation power and positively correlated with gamma oscillation power [14,16–18]. During movie watching, the fMRI signal in the human auditory cortex was found to be positively and negatively correlated to gamma and alpha band neuronal oscillations, respectively [19]. However, because this previous intracranial work utilized standard-resolution fMRI, it could not resolve how these dynamics are organized across cortical depths. Conversely, recent advances using high-resolution 7T fMRI have elucidated depth-specific neurovascular coupling in the visual cortex using simultaneous scalp EEG [20]. Yet, scalp-EEG suffers from volume conduction and signal attenuation, severely limiting its ability to resolve the highly localized HFA. Consequently, how precise intracranial neural dynamics and vascular responses are coupled across cortical depths in the human auditory cortex remains elusive. Here, we bridge this critical gap by combining the high spatial resolution of 7T fMRI (0.8 mm isotropic) with the spectral precision of stereoelectroencephalography (SEEG) during naturalistic music perception.

Previous animal studies have revealed that neuronal activation in the primary auditory cortex differs across cortical depths in response latency [21,22] and frequency selectivity [22,23]. The representation of auditory information along the vertical axis of the auditory cortex also differs across cortical layers: the granular layers have a more accurate tonotopic response than supragranular and infragranular layers [24–26]. A recent high-resolution fMRI study in the human auditory cortex suggest that the frequency preference is highly constant across cortical depths while task demands can sharpen the frequency tuning in the superficial depth more than the intermediate and deep depths [27]. However, how fMRI signals at different cortical depths are correlated to neural signals in the human auditory cortex is less clear.

Here, we use the temporal envelopes of oscillatory activity across distinct frequency bands, derived from intracranial neuronal recordings in participants under presurgical monitoring, to predict depth-specific hemodynamic responses to the same naturalistic stimulus in an independent dataset acquired using 7T fMRI [20,28,29] in healthy participants (see below for discussion of potential biases of sample characteristics). We aim to reveal neurovascular couplings across cortical depths to infer feedforward and feedback processing of auditory perception. With preferential anatomical feedforward connections to the intermediate cortical depth in the primary sensory region [30,31], we hypothesize that the hemodynamic response in intermediate cortical depth is more correlated to neural signals than in superficial or deep depths in processes where feedforward processing outweighs feedback. Specifically, alpha- and beta-band oscillations implicating feedback modulations are stronger in adaptation and higher cortical hierarchy, while neural activities at higher frequencies (> 30 Hz), including both oscillatory and non-oscillatory signals, relate to feedforward processing [32] of acoustic features.

Results

Figure 1 shows the locations of implanted contacts around the auditory cortex. In this study, we included measurements on all contacts within 15 mm of the centroid of the left and right auditory cortex defined by an atlas derived from the Human Connectome Project [33]. This range was suggested by our previous SEEG study, which showed the correlation among electrode contacts separated within 20 mm [34]. Together, six and seven patients had contacts within 15 mm to the centroid of left (n=33) and right (n=29) auditory cortex, respectively (Table S1). Three (2 patients) and three (2 patients) contacts were within 5 mm to the centroid left and right auditory cortices, respectively. None of the electrode contacts was close to the epileptic spike generators, as confirmed by the neurologists. The atlas also provided the delineation of primary (A1) and secondary (A2) auditory cortices (Figure 1). Because normalized depth (n.d.) is defined computationally via equidistant surfaces, these measurements represent geometric coordinates rather than absolute cytoarchitectonic boundaries. For clarity, we report our functional findings across ‘depths’, utilizing the BigBrain histological projection (Figure 1C) to provide qualitative anatomical context. Except n.d.=0.1, each normalized depth consisted of two to three cortical layers in the auditory cortex, A1, and A2. The depth n.d.-0.5 included the largest proportion for layer IV: 36%, 61%, and 30% for the auditory cortex, A1, and A2, respectively.

Figure 1.

fMRI signals entrained by musical stimuli vary across cortical depths

Following the experiment design of a previous study [19], our SEEG-fMRI correlation study used a cross-participant design: SEEG was measured from patients and fMRI was obtained from healthy participants. We first studied how the fMRI signals in the auditory cortex are entrained by acoustical features. Specifically, the depth-specific fMRI signals were correlated with the modeled fMRI signal generated by convolving the envelope of the presented music and a HRF [35]. Figure 2A shows the estimated coefficients in the linear regression. All coefficients were significant (FDR corrected p<0.001). Figure 2B shows the distribution of significant correlation across cortical depths. We found the most extensive and strongest correlation at n.d. = 0.1 and 0.3. Between cortical depths, correlations around n.d.=0.1 and 0.3 were stronger than pial surface and deep depths, while acoustic envelopes and fMRI signals at the pial surface are more significantly correlated than at deep depths (n.d.=0.5 to 0.9 and white matter boundary).

Figure 2

At the same cortical depth, there was no significant difference in the acoustic-fMRI correlation between A1 and A2 (Figure 2B). We also compared the acoustics-fMRI correlations between cortical depths within a region (Figure 2C). Among pair-wise comparisons across cortical depths, A1 at mid-to-superficial depths (n.d. = 0.3, which primarily maps to layer III in our histological projection; Figure 1C) demonstrated a significantly stronger acoustics envelope-fMRI signal correlation than the pial surface. At A2, this correlation was also significantly stronger at these mid-to-superficial depths (n.d. = 0.1 and 0.3) compared to deeper regions. The relationship between acoustic envelopes and fMRI signals across cortical depths were fitted to a cubic spline to create models for auditory cortex, A1, and A2. These models included the coupling between stimulus acoustics and fMRI signal and potential bias caused by vascular signals.

Correlations between fMRI and SEEG signals

Figure 3A and 3B show the between-subject[19] correlation between fMRI BOLD signals in the auditory cortex and the predicted responses by the convolution between band-limited SEEG activity patterns and a canonical HRF [35] during music hearing at different cortical depths. Significant negative correlations were found in alpha and beta bands at superficial (n.d.=0.1) and intermediate (n.d.=0.5) depths (FDR-corrected p<0.05). The envelope of gamma-band oscillation higher than 30 Hz convolved with HRF was positively correlated with the BOLD signal in all depths. Between cortical depths, fMRI signals at the geometric midpoint (n.d.=0.5) were less negatively and less positively correlated with alpha-beta and gamma band oscillations than at the most superficial depth (n.d.= 0.1), respectively. Related to deep depth (n.d.=0.9), intermediate depth (n.d.=0.5) fMRI signals were more negatively and positively correlated with alpha and gamma band oscillations, respectively.

Figure 3

Distributions of significant correlations between depth-specific fMRI signals and neural oscillations in alpha, beta, and gamma bands are shown in Figure 3C. Much stronger and diffusive patterns in the temporal lobes were observed in the gamma band than in the alpha and beta bands, with a trend of increasing correlations toward superficial depths. The correlations in alpha and beta bands are prominent in the left hemisphere.

The correlation between fMRI signals and oscillatory SEEG signals in A1 and A2 was similar across cortical depths (Figure 3D). Both regions exhibited negative correlations in the alpha and beta bands and positive correlations in the gamma band. The statistical comparison (Figure 3D, bottom panel) revealed that A2 had significantly more negative correlations to alpha and beta band oscillations compared to A1. This difference was most prominent across deeper cortical depths (n.d.= 0.7 to 0.9), with a localized, secondary effect also observed at the pial surface. In the gamma band between 50 Hz and 150 Hz, the fMRI-SEEG signal correlation at A2 was significantly stronger than A1 at deep depths (n.d.=0.7, 0.9, and white matter boundary). A significantly more negative correlation between fMRI signals and low-frequency oscillations (4-6 Hz) at A1 compared to A2 was observed across most cortical depths, extending from the pial surface down to the white matter boundary.

Neural oscillations derived from more electrodes were distributed over a wider range (15 mm from the auditory cortex centroid). Therefore, the specificity of the neural oscillation estimates may be compromised in improving the signal by averaging data from more electrodes. We studied this effect by limiting electrode contacts (n = 6; three in the left hemisphere and three in the right hemisphere) only within 5 mm of the auditory cortex. Repeating the correlation between fMRI and the more selective SEEG measurements shows similar results (Supplementary Figure 1). Namely, stronger negative correlations were found at the alpha and beta bands and positive correlations at the gamma band. Again, the correlations became more significant toward superficial depths. Distributions of the positively correlated gamma oscillations and fMRI signals were more focal when we limited to electrode contacts within 5 mm of the auditory cortex centroid.

HFA correlates fMRI signals across cortical depths

HFA taken from electrode contacts within 5 mm of A1 had significant correlation to depth-specific fMRI with a trend of stronger correlation toward the pial surface (Figure 4A). The strongest correlation between fMRI signals within A1 was at n.d.=0.3, while A2 and the whole auditory cortex had stronger correlation at n.d. = 0.3, 0.1, and pial surface (Figure 4B). Contrasting the HFA-fMRI correlation with acoustics-fMRI correlation (Figure 2), the HFA-fMRI correlation in the auditory cortex was significantly stronger at intermediate-deep (n.d.=0.5, 0.7, and 0.9) but weaker at superficial (n.d. = 0.1) depths (Figure 4B). Specifically, A1 showed a selectively stronger HFA-fMRI correlation compared to its acoustic-fMRI correlation at the white matter boundary, the geometric midpoint (n.d.=0.5), and the superficial depth (n.d.=0.1). Between cortical depths, the most significant correlation difference was found between intermediate (n.d. =0.5) and deep depths (n.d.=0.9 and white/gray matter boundary) (Figure 4C). The centroid of the area showing significant HFA-fMRI signal correlation was found moving away from A1 to A2 toward the latero-posterior direction at deeper cortical depths (Figure 4D), suggesting superficial feedback type signaling increased across the auditory hierarchy from posteriolateral to anteromedial aspects of auditory cortex.

Figure 4

Relationships of neuronal oscillations and fMRI signals across cortical depths change between repeated listening

As we measured SEEG and fMRI with repeated music presentation, we examined how the relationship between oscillatory neural and fMRI signals across cortical depths changes between the first and second listening. Supplementary Figure 2 shows these relationships at both A1 and A2. Visually, contrasting between the first and second listening, we found stronger negative correlations between alpha/beta band oscillations and fMRI signals and stronger positive correlations between neural activity above 70 Hz up to 200 Hz and fMRI signals in both A1 and A2. The significant differences were found more across frequencies, particularly between 50 Hz and 200 Hz, in A2 than A1. At A1, more significantly negative correlations in alpha and beta bands were found at intermediate (n.d.=0.5) depths. At A2, more significantly negative alpha/band-fMRI correlations were observed across cortical depths. Significantly stronger correlations at broadband high frequency activity (> 70 Hz) and fMRI signals were found across cortical depths with stronger effects above n.d.=0.5 at A1 selectively, while this difference was found across all depths at A2. Interestingly, we also found that the correlation between low gamma band (between 30 Hz and 50 Hz) oscillations and fMRI were reduced in the repeated listening, particularly around intermediate and superficial depths, at both A1 and A2. At the same time, a stronger correlation between theta band (between 4 Hz and 8 Hz) oscillation and fMRI signals was found at A1 and A2.

The consistency between within-subject and between-subject fMRI-SEEG correlation

To study whether the correlation between fMRI and SEEG signals from the same group of subject is different from that using SEEG and fMRI signals from different groups, we averaged all 7T data across cortical depths from healthy controls and calculated the frequency-dependent fMRI-SEEG correlation. The result was compared with the epilepsy patients fMRI data from 3T. Supplementary Figure 3 shows agreement between correlations using data at 7T and 3T in intra-subject and inter-subject analysis, particularly the positive correlation between gamma-band oscillations and fMRI signals. These results provide evidence supporting the validity of the correlation analysis across participant groups.

Discussion

We studied the relationship between intracranial neuronal activity and depth-specific fMRI signals in the human auditory cortex during perception of complex naturalistic stimuli. Leveraging invasive recordings from epilepsy patients, we achieved high spatial resolution measurements of local field potentials. Combined with ultrahigh-resolution 7T fMRI (0.8 mm isotropic), this approach enabled us to characterize neurovascular coupling across primary and secondary auditory cortices during music processing.

The correlation between acoustic features and fMRI signals spanned all cortical depths, consistent with recent laminar recordings in the mouse auditory cortex [36]. The strongest absolute correlations between fMRI signals and the envelope of musical stimuli were found at mid-to-superficial depths (n.d.=0.3 and 0.1) (Figure 2), likely reflecting the general vascular draining bias toward the cortical surface. Because this macroscopic vascular bias obscures the underlying neural directionality, we leveraged our intracranial recordings to investigate how distinct, frequency-specific neural dynamics drive these depth-dependent hemodynamic responses.

In neural oscillation analysis, we found that fMRI signals correlate negatively with alpha/beta-bands and positively with the gamma band. These correlations are more prominent at superficial depths than at the deep white-gray matter boundary (Figure 3). This is consistent with a correlational study between single neuron recordings and fMRI during music perception [19]. The coupling between band-limited neural oscillations in alpha, beta, and gamma bands and fMRI signals did not show selective peaks at intermediate, superficial, or deep depths. Instead, a strong correlation (positive in gamma and negative in alpha/beta bands) was found at the superficial depth with a diminishing gradient toward deeper depths. Together with the selective enhancement of HFA-fMRI correlations relative to acoustics-fMRI correlations at intermediate depths (Figure 4), our frequency-specific results demonstrate that naturalistic music perception seamlessly integrates both feedforward and feedback processing streams. This spectral dissociation aligns with predictive coding frameworks in auditory processing, where gamma activity is associated with the feedforward propagation of sensory inputs or prediction errors, while beta oscillations mediate top-down feedback predictions [37].

The neurophysiological basis of fMRI signal has been attributed to local field potentials in non-human primate studies [13,14]. These neural activities coupled to fMRI signals were represented in oscillations faster than 50 Hz [14,19]. Notably, HFA (>70 Hz) has traditionally been taken as a proxy for asynchronous multi-unit activity and local neuronal spiking [38,39]. However, recent evidence suggests that HFA is not a monolithic signal; it comprises distinct early components related to local spiking and late components driven by dendritic processes [7]. Our HFA results (Figure 4) highlight the complexity of this relationship across the cortical column. Anatomically, ascending feedforward pathways predominantly target intermediate layer IV. The robust HFA-fMRI correlation we observe peaking at these intermediate depths likely reflects the strong, early spiking activity driven by bottom-up acoustic inputs. Conversely, the spread of significant HFA-fMRI correlations to superficial layers —which are heavily targeted by top-down cortico-cortical connections [40,41] —suggests that depth-specific fMRI also captures the late, dendritic-driven components of HFA associated with feedback modulation. Therefore, rather than segregating spiking and synaptic activity by depth, our data suggests that the intermediate-depth HFA correlation reflects dominant feedforward sensory driving, while superficial correlations capture subsequent modulatory processing. Note that the link between synaptic input and fMRI signals was also supported by the anti-correlation in alpha and beta bands [2] (Figure 3). Taken together, our results using cortical depth-specific fMRI and invasive neural recording suggest a tight coupling between feedforward and feedback modulations involving both neural input and output activities during music perception.

As we move up the cortical hierarchy, feedback modulation plays an increasingly prominent role in information processing. To investigate this, we examined how the coupling between band-limited neural oscillations and cortical depths varies across the auditory hierarchy human auditory cortex by comparing A1 and A2 (Figure 3). We found that A2 has a stronger coupling between fMRI signals and alpha and beta oscillations at superficial (n.d.=0.1) and deep (n.d.=0.7 and 0.9) depths than A1. Considering the linkage between alpha/beta oscillations and feedback modulation [42,43] and the feedback modulations at infragranular and supragranular layers [2,3,44], our results support that more feedback modulation is present at locations of a higher cortical hierarchy in the auditory perception.

The neurovascular coupling studied here exhibited a general trend of a stronger correlation at superficial depths, even with our efforts of regressing out the average fMRI signal at pial and gray/white matter boundaries. This can be attributed to the venous flow from deep to superficial depths [45]. This feature has also been observed in EEG-fMRI correlations in the human visual cortex [20]. However, some of our analyses also suggest depth-specific properties: in pair-wise comparisons of fMRI-acoustics (Figure 2C) and fMRI-SEEG HFA correlation (Figure 4C), the largest difference was found between intermediate-superficial depths (n.d.=0.5, 0.3, and 0.1) and white matter boundary. Furthermore, in the comparison between repeated listening (Figure S1), venous bias was also reduced because this within-region vascular bias was not likely to change between conditions[46]. The reduced negative coupling between alpha/beta oscillations and fMRI signals may reflect a mechanism that prioritizes auditory processing [47]. The increased coupling between theta oscillations and fMRI signals may follow from reduced cognitive load during repeated listening [48]. Similar strategy of contrasting between conditions was used in taking fMRI-acoustic correlation as a reference: the fMRI-SEEG HFA correlation differs from fMRI-acoustics correlation mostly at deep and superficial depths (Figure 4B). More complicated methods, such as biophysical modeling [49], may suppress the venous bias more specifically.

Depth-specific fMRI has been used most frequently in visual [50–52] and somatosensory cortices [53]. In auditory cortices, it has been shown that the tonotopic mapping remains stable across cortical depths, while attention can tune acoustic frequency preference more than intermediate and deep depths [27]. Our results show how depth-specific fMRI signals are correlated with alpha, beta, and gamma band oscillations across the auditory cortex hierarchy. Together, we showcase how to use sub-millimeter fMRI to enrich the understanding of the neurophysiological underpinning of auditory perception and cognition.

A limitation of our approach is the translation of computationally defined normalized cortical depths to true cytoarchitecturally-defined cortical layers. While we qualitatively matched our equidistant surfaces to the BigBrain dataset [54] —showing, for instance, that layer IV is most heavily represented at n.d.=0.5 —cortical thickness and curvature variations mean that a depth like n.d. = 0.3 captures a blend of supragranular and granular features depending on the specific location (Figure 1C). Future quantitative studies using subject-specific high-resolution quantitative MRI are needed to achieve exact histological correspondence.

This study used musical stimuli, a potent acoustic stimulus known to influence emotions through both the auditory and limbic systems [55], to investigate how neuronal activity and blood flow are linked across different depths of the auditory cortex. This design echoes the increasing interested in exploring human brain function with more ecologically relevant stimuli [56]. Our stimuli included instrumental and vocal components in different genres, leading to the cortical-depth specific correlations that are expected to be generalized to other musical stimuli.

Our study used SEEG and fMRI data from different participants. This between-subject design [19] carries the advantage of avoiding the carryover effect (between SEEG and fMRI) but cannot control the tentative confound attributed to the difference between healthy participants providing fMRI measurements and epilepsy patients providing SEEG recordings. Taking SEEG and fMRI from different individuals was due to the practical challenges of scheduling epilepsy patients for 7T MRI. Note that to reduce the concern of effect from groups, all SEEG data from contacts showing epileptic features were excluded by clinicians. Future studies with SEEG and fMRI on the same patient can further address this concern.

Taken together, our findings demonstrate that depth-specific fMRI, when combined with intracranial electrophysiology, can reveal a spectrolaminar organization of neurovascular coupling in the human auditory cortex during naturalistic perception. By linking frequency-specific neuronal dynamics to hemodynamic responses across cortical depths, we provide novel insights into the laminar architecture of feedforward and feedback processing in humans. This approach not only bridges invasive and non-invasive modalities but also opens new avenues for investigating the neural basis of complex cognitive functions and their disruption in neuropsychiatric disorders.

Materials and Methods

All epilepsy patients and healthy participants joined this study with written informed consent after the approval of the Institute Review Board from Taipei Veteran General Hospital and SungKyunKwan University, respectively. All participants had normal hearing. We presented three songs (song 1: Japanese animation Doraemon theme, 176 s; song 2: Brahms Piano Concerto No. 1, 360 s; Song 3: Lost Stars from Adam Levine, 266 s) to each individual during fMRI or stereoelectroencephalography (SEEG) recordings using Psychtoolbox in Matlab. These musical stimuli are available at https://github.com/fahsuanlin/labmanual/wiki/32:-Sample-data:-SEEG-recording-during-music-listening. Each participant was instructed to relax and listen to musical pieces, each of which was played twice in a randomized order. A 10-s silence interval was inserted as rest condition between two songs, before the first song, and after the last song.

For healthy participants, functional MRI data were acquired on 7T MRI (Terra, Siemens Healthineers, Erlangen, Germany) with a 32-channel whole-head coil array (Nova Medical, Wilmington, MA, USA). Structural and functional images, both with 0.8 mm isotropic resolution, were acquired with MP2RAGE [57] and gradient-echo EPI (TR=2.8 s, TE=29 ms, flip angle=75°, 64 slices), respectively. The fMRI covered the temporal lobe at both hemispheres. The fMRI acquisition was “blip-down”, where the phase encoding direction was from the positive to the negative direction. For each participant, we additionally acquired EPI “blip-up” images, which had phase encoding from negative to positive phase encoding direction. The combination of “blip-down” and “blip-up” trajectories allowed for estimating the field map, which was used to correct image geometric distortion (see description below). Twelve epilepsy participants had MRI data acquired on 3T MRI (Skyra, Siemens Healthineers, Erlangen, Germany) with a 32-channel whole-head coil array. Structural images in 1-mm isotropic resolution were acquired with MPRAGE sequence and functional images in 1.5-mm isotropic resolution were acquired with gradient-echo EPI (TR=2.5 s, TE=30 ms, flip angle=60°, 38 slices).

Structural images of healthy participants were reconstructed to delineate the gray-white matter boundary automatic segmentation and cortical parcellation using FreeSurfer ^9,10 (http://surfer.nmr.mgh.harvard.edu; version 7.2). Nine cortical surfaces with equally spaced cortical thickness were reconstructed (mris_expand in FreeSurfer). In this study, we used “normalized depth” (n. d.) to denote the interpolated surface. For example, n. d. = 0.1 was the surface about 10% to the gray-white matter boundary and about 90% from the pial surface. The relationship between normalized depths to cortical layers was quantified by morphing cortical layer thickness [54] from the BigBrain data, a 3D histological model of the human brain with 20 mm resolution [58], to FreeSurfer atlas using the BigBrainWarptoolbox[59] (https://github.com/caseypaquola/BigBrainWarp). While the equi-volume principle is theoretically advantageous for accounting for cortical folding in cytoarchitecture, we opted for the equi-distance approach (mris_expand). At our in-vivo spatial resolution (0.8 mm isotropic), the equi-volume algorithm can severely amplify discretization and segmentation noise due to its sensitivity to local curvature estimates. Furthermore, empirical comparisons at this resolution demonstrate that partial voluming effects dominate, rendering the differences between equi-volume and equi-distance layering negligible for functional outcomes [60].

Functional MRI data from healthy participants were first preprocessed by motion correction and slice timing correction using the FS-FAST pipeline in FreeSurfer. Image distortion due to off-resonance was corrected by estimating the field maps from “blipped-up” and “blipped-down” images[61] (topup function in FSL; v. 5.0; https://fsl.fmrib.ox.ac.uk/fsl/docs/). Distortion-corrected fMRI was spatially registered to the structure image using boundary-based registration [62] (bbregister in FreeSurfer). Finally, fMRI data were estimated on each depth-specific surface using mri_vol2surf function in FreeSurfer. fMRI data were estimated over the pial surface (gray matter and cerebrospinal fluid boundary), white/gray matter boundary (white surface in this study), and nine equispaced intermediate surfaces between (n.d.=0.1 to n.d=0.9 in steps of n.d.=0.1). No spatial smoothing was applied between depths. The specificity of cortical depth fMRI is degraded by venous signals in gradient-echo EPI in amplifying the signal toward the cortical surface [63]. To reduce this bias, we first calculated the average fMRI signals at all surfaces between the pial and white/gray matter boundary. Then, a linear regression was used on the depth-specific fMRI signal to remove this average signal at each cortical location.

Electrophysiological responses were measured from 16 medically refractory epilepsy patients by SEEG. Locations of the implanted electrodes were planned based on each patient’s clinical need. All patients had an electrode (Ad-Tech Medical Instrumnt, Oak Creek, WI, USA) with up to ten contacts (5 mm separation) at the temporal lobe. All electrode contacts showing epileptic-like activity verified by neurologists were excluded from the analysis. We excluded data from three patients because they had no contacts within 15 mm from the centroid of the left or right auditory cortex defined from an atlas [33]. Eventually, SEEG data were utilized for further analysis from 13 patients and 62 electrode contacts.

Pre-surgery and post-surgery structural MRI’s (3D T₁-weighted images with MPRAGE sequence; 3T MRI, Skyra, Siemens Healthineers, Erlangen, Germany; 1 mm isotropic resolution) were obtained from patients to identify electrode and contact locations [64].

The SEEG data were re-referenced to an electrode contact at the white matter away from any tentative region showing epileptic activity. Frequency-specific oscillatory neural activities were estimated by first applying the Morlet wavelet transform (the central frequencies varying between 4 Hz and 300 Hz and 7-cycle width; 1 Hz step below 10 Hz, 2 Hz step below 20 Hz, 5 Hz step below 60 Hz, 10 Hz step below 120 Hz, 20 Hz step above 120 Hz) to the SEEG time series and then taking the absolute values. At each central frequency, a modeled fMRI time series was created by convolving a canonical hemodynamic response function (HRF) [35] to the oscillatory neural activities at that frequency. The General Linear Model (GLM) was used to correlate between frequency-specific oscillatory neural activities and depth-specific fMRI data. Specifically, in the band-specific calculation, all SEEG-modeled fMRI time series at frequencies within alpha (8 Hz – 14 Hz), beta (14 Hz – 30 Hz), low gamma (30 Hz – 60 Hz), and high gamma (60 Hz – 300 Hz) bands were included as GLM regressors. These neural oscillations characterized the local field potentials.

Broadband high-frequency activity (HFA) spanning between 70 Hz and 150 Hz reflects non-rhythmic asynchronous neural firing linked to local cortical processing [38,39]. We also examined how HFA is correlated with depth-specific fMRI signals: The envelop of band-pass (70 – 150 Hz) filtered SEEG signals was extracted by taking absolute values after Hilbert transform. These HFA from electrode contacts at the vicinity of the auditory cortex were first averaged and then convolved with a canonical HRF [35]. This convolved signal was taken as a regressor in GLM to correlate between HFA and depth-specific fMRI data.

To study how acoustic features correlate with fMRI data, we also extracted the envelope of the presented music using Hilbert transform (envelope function in Matlab, Natick, MA, USA). The extracted envelope was convolved with the same canonical HRF to model the fMRI time series, which was used to correlate with empirical depth-specific fMRI data using GLM. We used cubic spline to build a model describing how the correlation between fMRI signals and acoustic envelope changes across cortical depths. This model was fitted to the HFA-fMRI correlation to contrast the difference between neural and acoustics correlation to fMRI across cortical depths.

In all GLM, we included six motion parameters and time series in the ventricles and white matter segmented by FreeSurfer as confounds. We also include a constant and linear trend to model the nuisance effects. We first obtained the regression coefficients from each participant. Subsequently, we performed a second-level GLM to evaluate the signficance across participants across songs and repeated listening.

We particularly focused on the primary (A1) and secondary (A2) auditory cortex defined by the Human Connectome Project data [33]. All statistical tests were subjected to controlling the false discovery rate (FDR) to avoid the inflation of type-I error in multiple comparisons.

To study the consistency between correlating within-subject and between-subject, we correlated epilepsy patients’ 3T fMRI with SEEG oscillatory power envelopes across frequencies. At the same time, we also averaged healthy controls’ 7T data across cortical depths and correlated with the same oscillatory power envelopes of SEEG.

Supplementary material

Fig S1

Fig. S2

Fig. S3.

Table S1.

Data availability

Data Availability Due to ethical and privacy restrictions governing human patient data-specifically, intracranial stereoelectroencephalography (SEEG) recordings from epilepsy patients and high-resolution 7T MRI scans-the raw data analyzed in this study cannot be made publicly available. De-identified data can be shared with qualified researchers for non-commercial research purposes upon formal request to the corresponding authors and the Institutional Review Board (IRB) of Taipei Veterans General Hospital. Access will be granted subject to ethical approval and the establishment of a formal Data Transfer Agreement. Code Availability All custom MATLAB scripts and analytical pipelines used to process the fMRI and SEEG data, including the layer-specific extraction and cross-correlation analyses, are publicly available on GitHub at https://github.com/fahsuanlin/labmanual/wiki.

Acknowledgements

We thank funding from the following institutes to support this research: Canadian Institutes of Health Research grant PJT 178345 and PJT 496433 (FHL), National Sciences and Engineering Research Council grant RGPIN-2020-05927 (FHL), Canada Foundation for Innovation grant 38913 and 41351 (FHL), MITACS IT25405 (FHL).

Additional information

Author Contributions

Conceptualization: HYY, WJK, FHL

Methodology: HJL, HYY, KU, WJK, FHL

Investigation: HJL, JA, FHL

Visualization: FHL

Supervision: WJK, FHL

Writing—original draft: HJL, JA, FHL

Writing—review & editing: HJL, JA, HYY, CCC, CCL, WJK, HL, KU, FHL

Funding

Canadian Institutes of Health Research (CIHR) (PJT 178345)

• Fa-Hsuan Lin

Canadian Institutes of Health Research (CIHR) (PJT 496433)

• Fa-Hsuan Lin

Natural Sciences and Engineering Research Council of Canada (NSERC) (RGPIN-2020-05927)

• Fa-Hsuan Lin

Mitacs (Mitacs Inc.) (IT25405)

• Fa-Hsuan Lin

Canada Foundation for Innovation (CFI) (38913)

• Fa-Hsuan Lin

Canada Foundation for Innovation (CFI) (41351)

• Fa-Hsuan Lin