Static and Dynamic Resting‐State fMRI Connectivity Show Frequency‐Resolved Associations With Intracranial EEG Connectivity

2.1. Participants

Seventy individuals undergoing intracranial video‐EEG monitoring (VEM) for drug‐resistant focal epilepsy were prospectively recruited between 2008 and 2021 for simultaneous icEEG‐fMRI acquisition. Inclusion criteria were age > 18 years, ability to provide informed consent, absence of MRI contraindications, and a clinical diagnosis of focal epilepsy. The study was approved by the University of Calgary Conjoint Health Research Ethics Board.

2.2. Intracranial EEG Data Acquisition

Each subject's intracranial VEM data were reviewed by an experienced epileptologist (P.F.) to select up to eight electrodes to be used during the icEEG‐fMRI acquisition. Selected electrodes were connected through an MRI‐compatible electrode connector system supporting up to 64 channels (Compumedics Neuroscan) and routed to a SynAmps RT headbox and amplifier system (Compumedics Neuroscan). EEG signals were acquired at 10 or 20 kHz using Scan 4.4 software before 2020 or Curry 8 software from 2020 onward. The EEG system clock was synchronized with the MRI scanner clock. Online artifact rejection and filtering were used to facilitate EEG monitoring during acquisition.

2.3. MRI Data Acquisition

Subjects were scanned using a 3 T GE Discovery MR750 whole body scanner with an eight‐channel receive‐only phased‐array head coil (GE Healthcare). Subjects were closely monitored during the acquisition. The MRI protocol included multislice anatomical imaging [spoiled gradient‐recalled echo two‐dimensional (2D) multislice sequence: echo time = 2.1 ms, repetition time = 150 ms, flip angle = 18°, 128 × 128 matrix, 24 slices 5.00 mm thick], anatomical 3D T1‐weighted imaging (echo time = 3.8 ms, repetition time = 9.3 ms, flip angle = 12°, 24‐cm field of view, 320 × 256 × 64 matrix, 0.47 × 0.47 × 2.00 mm slices), and fMRI (gradient‐recalled echo planar imaging: echo time = 30 ms, repetition time = 1500 ms, flip angle = 65°, 24‐cm field of view, 64 × 64 matrix, 24, 25, or 29 slices 3.75 × 3.75 × 5.00 mm). During resting‐state fMRI acquisition, participants were instructed to keep their eyes closed. Subjects were scanned for up to 60 min, comprising three 20‐min resting‐state fMRI runs.

2.4. Intracranial EEG Preprocessing

icEEG data were preprocessed using an in‐house MATLAB pipeline (R2022b; MathWorks, Natick, MA). For the purpose of gradient artifact correction, channels were re‐referenced using a within‐electrode common average to minimize gradient artifact amplitude while preserving the flexibility of a referential montage for subsequent identification of epileptiform activity. The preprocessing pipeline included: (i) bandpass filtering (0.5–500 Hz); (ii) correction of MRI gradient artifacts using average artifact subtraction (Allen et al. 2000); (iii) reduction of residual artifacts using principal component analysis across TR epochs (Niazy et al. 2005); and (iv) removal of line noise (60 Hz and harmonics) using notch filtering. Data were subsequently downsampled to 2500 Hz.

2.5. fMRI Preprocessing

All MRI data were processed using FEAT version 6.0, part of the FSL FMRIB Software Library (FSL; https://fsl.fmrib.ox.ac.uk/fsl/docs/). Standard preprocessing included slice‐timing correction, brain extraction (Smith 2002), head motion correction (Jenkinson et al. 2002), spatial smoothing with a 6‐mm FWHM Gaussian kernel, and high‐pass temporal filtering with sigma = 50 s. Frames with head motion exceeding 1.5 mm were removed, and affected runs were split into contiguous segments at censored frames. Segments were not concatenated across motion boundaries; rather each retained segment was treated as a separate run. Segments shorter than 5 min were excluded from further analysis.

FSL's multivariate exploratory linear optimized decomposition into independent components (MELODIC) tool was used to decompose each functional scan into 60 components. All components were visually inspected by two reviewers (T.R., W.W.) to identify noise‐ and or artifact‐related components (Griffanti et al. 2017), which were then regressed out of the data. Functional MRI data were registered to both low‐ and high‐resolution anatomical images using FMRIB linear image registration tool (FLIRT) (Jenkinson et al. 2002; Jenkinson 2001).

2.6. Regions of Interest (ROI) and Electrode‐to‐fMRI Mapping

Electrode contacts were manually localized on each subject's post‐implantation T1‐weighted image. Bipolar channels were defined as differences between adjacent contacts along each depth electrode and served as the basis for all analyses. Channels exhibiting interictal epileptiform discharges (IEDs) were excluded. For each remaining bipolar channel, a corresponding fMRI ROI was defined in native space by combining 3 × 3 × 3 voxel neighborhoods centered on the two constituent contacts, from which the mean BOLD time series was extracted and bandpass‐filtered between 0.01 and 0.1 Hz. Pairwise FC was then computed between bipolar icEEG channels and between their corresponding fMRI ROIs. Throughout the manuscript, a sampled electrode‐pair connection refers to a pair of bipolar icEEG channels together with the corresponding pair of fMRI ROIs. To prevent inflated correlations due to overlapping anatomical sampling, connections between bipolar channels sharing a common contact were excluded from all coupling analyses.

2.7. Functional Connectivity Analyses (icEEG)

Static and dynamic icEEG FC were computed to characterize time‐averaged and window‐level interactions during rest. Continuous icEEG signals were filtered into canonical frequency bands: delta: 1–3 Hz; theta: 4–7 Hz; alpha: 8–12 Hz; beta: 13–30 Hz; low gamma: 31–50 Hz; high gamma: 51–100 Hz. While higher frequencies are available in icEEG, the upper limit of 100 Hz was chosen a priori to balance sensitivity to high‐frequency icEEG activity with signal quality and interpretability in simultaneous icEEG–fMRI recordings, as higher‐frequency activity is more susceptible to lower signal‐to‐noise ratio due to residual gradient artifacts and could also be confounded by higher‐frequency oscillations (HFO)–range activity (Frauscher and Gotman 2019). This upper limit was selected a priori to balance sensitivity to high‐frequency icEEG activity with signal quality and interpretability in simultaneous icEEG–fMRI recordings. Band‐limited power (BLP) time series were obtained using the Hilbert transform in MATLAB. Full‐recording static icEEG FC was defined as the Pearson correlation between BLP time series for each pair of bipolar channels across the full retained recording. Dynamic icEEG FC was estimated using 120 s sliding windows (Leonardi and Van De Ville 2015) with 60 s overlap (18 windows per run), with windowed correlations computed for each channel pair to obtain time‐resolved FC estimates (Figure 1). A window‐averaged icEEG FC was also computed by averaging sliding‐window FC estimates across retained windows and runs, providing a whole‐scan summary based on the same windowed estimator used in the dynamic analyses.

2.8. Functional Connectivity Analyses (fMRI)

Static and dynamic resting‐state fMRI FC were computed using ROIs corresponding to bipolar icEEG locations. BOLD time series were extracted from ROIs corresponding to bipolar icEEG locations and shifted by 6 s, corresponding to 4 TRs, to approximate the expected hemodynamic delay. Full‐recording static fMRI FC was defined as the Pearson correlation between ROI BOLD time series across the full retained recording. Dynamic fMRI FC was estimated using the same sliding‐window framework as the icEEG analysis, with 120 s windows and 60 s overlap (Figure 2). A window‐averaged fMRI FC summary was also computed by averaging sliding‐window fMRI FC estimates across retained windows and runs. Correlation values were Fisher z‐transformed before statistical analysis.

2.9. Static icEEG–fMRI Coupling Analysis

To assess stable correspondence between electrophysiological and hemodynamic connectivity, static icEEG–fMRI coupling was evaluated using the two FC summaries defined above. First, full‐recording static FC was computed over the full retained recording duration for each subject, electrode‐pair connection, and frequency band, yielding a conventional whole‐recording measure of time‐averaged connectivity. Second, window‐averaged FC was computed as the mean of sliding‐window FC estimates across retained windows and runs, yielding one mean connectivity value per modality for each electrode pair. Although derived from a windowed framework, this estimator represents a static summary of connectivity and was included to provide direct comparability with the dynamic window‐level analyses.

2.10. Dynamic Window‐Level icEEG–fMRI Coupling Analysis

To assess window‐level coupling between electrophysiological and hemodynamic connectivity, sliding‐window icEEG FC was decomposed into two components: (i) a subject‐, electrode‐pair‐, and frequency band‐specific mean icEEG FC term, computed as the average of sliding‐window icEEG FC estimates across retained windows and runs; and (ii) a window‐specific icEEG FC deviation term ($\text{ΔicEEG}$), computed as the difference between each windowed icEEG FC estimate and the corresponding mean sliding‐window icEEG FC. In other words:

This decomposition was used to separate average differences in icEEG FC across sampled electrode‐pair connections from fluctuations around that average within the same connection. Thus, ΔicEEG was centered by construction within each subject, electrode‐pair connection, and frequency band, and was interpreted as a measure of within‐connection window‐level variation in icEEG FC.

Specifically, window‐level icEEG FC was expressed as follows:

where ${\overline{\text{icEEG}\_\mathrm{FC}}}_{s,e,b}^{\mathrm{win}}$ denotes the mean sliding‐window icEEG FC for subject $s$, electrode‐pair connection $e$, and frequency band $b$, computed across retained windows and runs. $\mathrm{\Delta }{\text{icEEG}\_\mathrm{FC}}_{s,e,b,r,w}^{\mathrm{win}}$ denotes the window‐specific deviation from that mean for subject $s$, electrode‐pair connection $e$, frequency band $b$, run $r$, and window $w$. Windowed fMRI FC for the corresponding electrode‐pair connection and time window was then modeled as a function of both the mean icEEG FC term and the window‐specific ΔicEEG term, allowing average electrophysiological coupling across sampled connections to be distinguished from window‐level electrophysiological fluctuations. Because the mean icEEG FC term and ΔicEEG were derived from the same sliding‐window icEEG FC estimates, they were not interpreted as independent measurements. Rather, this decomposition was used to partition the same windowed icEEG FC signal into an average connection‐specific component and a window‐specific deviation component. Full model equations and random‐effects specifications are provided in the Supporting Information.

2.11. Linear Mixed‐Effects Modeling and Statistical Inference

Associations between electrophysiological and hemodynamic connectivity were tested using linear mixed‐effects models (LMMs). Models were used to account for the hierarchical structure of the data, including repeated observations within subjects and across electrode‐pair connections, and for dynamic analyses, repeated observations across time windows.

For static icEEG–fMRI coupling, models were fitted separately for each frequency band to test whether variation in icEEG FC across sampled electrode‐pair connections was associated with corresponding variation in fMRI FC. Separate models were fitted for the full‐recording static FC estimator and the window‐averaged FC estimator. Fixed effects included the relevant icEEG FC term, icEEG FC variability where applicable, inter‐electrode distance, age, and sex. Inter‐electrode distance was log‐transformed and mean‐centered to account for systematic distance‐related effects on connectivity estimates. Random effects included subject‐ and electrode‐pair‐level terms to account for repeated observations and inter‐individual variability.

Because the dynamic analyses were performed at the window level, we first characterized the variance structure of windowed FC using intraclass correlation coefficients (ICCs). ICCs were computed separately for windowed fMRI FC and band‐limited windowed icEEG FC using random‐intercept models with subject as the grouping factor. For fMRI, duplicate rows across icEEG frequency bands were removed because the same windowed fMRI FC estimates were repeated for each band. ICCs were interpreted as the proportion of total window‐level FC variance attributable to between‐subject differences. Bootstrap standard deviations were estimated by resampling subjects with replacement 500 times.

For dynamic window‐level coupling, band‐wise models were fitted separately for each frequency band, with windowed fMRI FC for each electrode‐pair connection and temporally aligned time window as the dependent variable. Fixed effects included the mean sliding‐window icEEG FC term, the window‐specific icEEG FC deviation term (ΔicEEG), and icEEG FC variability, quantified as the standard deviation of windowed icEEG FC across retained windows and runs for each subject, electrode‐pair connection, and frequency band. Inter‐electrode distance, age, and sex were included as covariates. Random effects included electrode‐pair connection, run nested within subject, and subject‐specific random slopes for ΔicEEG. Statistical inference focused on the fixed effect of $\text{ΔicEEG}$, which tested whether within‐connection, window‐level deviations in icEEG FC were associated with temporally aligned windowed fMRI FC after accounting for mean icEEG FC, overall icEEG FC variability, inter‐electrode distance, age, and sex.

To formally test whether the $\text{ΔicEEG}$ coupling effect differed across frequency bands, we fitted an additional omnibus dynamic LMM including frequency band and a band × $\text{ΔicEEG}$ interaction term. This model was compared with a reduced model without the interaction term using a likelihood‐ratio test. Models were fitted using maximum likelihood for this comparison. This analysis tested whether the fixed effect of $\text{ΔicEEG}$ differed across frequency bands, rather than relying only on differences in FDR‐corrected significance across separate band‐wise models.

Multiple comparisons across icEEG frequency bands were controlled using the Benjamini–Hochberg false discovery rate procedure within each analysis. Reported coefficients reflect fixed‐effect estimates from the corresponding LMMs. Full model equations, random‐effects specifications, ICC analyses, and sensitivity/control analyses are provided in the Supporting Information. An overview of the analysis is shown in Figure 3.

2.12. Sensitivity and Robustness Analyses

Several sensitivity and robustness analyses were performed to assess whether the reported icEEG–fMRI associations were robust to key preprocessing choices, dynamic FC parameters, anatomical sampling, and head motion. First, to assess whether the results were influenced by the 6‐mm FWHM spatial smoothing kernel used in the primary fMRI preprocessing pipeline, we repeated the icEEG–fMRI coupling analyses using fMRI data processed with alternative 4‐mm and 8‐mm FWHM spatial smoothing kernels. All other preprocessing and analysis parameters were kept identical to the primary analysis. For the static analyses, we assessed whether full‐recording static FC and window‐averaged FC associations were robust across smoothing kernels. For the dynamic analyses, statistical inference focused on the primary window‐level coupling term, ΔicEEG, while mean sliding‐window icEEG FC and icEEG FC variability were retained as covariates in the model. The 4‐mm and 8‐mm dynamic analyses were restricted to the same subject, run, sampled electrode‐pair connection, frequency band, and window keys used in the primary 6‐mm analysis. These analyses were used to assess whether the reported icEEG–fMRI coupling results were sensitive to the choice of fMRI spatial smoothing kernel.

Second, to assess sensitivity to temporal alignment between icEEG and fMRI connectivity estimates, we repeated the dynamic icEEG–fMRI coupling models using fMRI temporal shifts of 0, 3, and 9 s, in addition to the primary 6 s shift. The 6 s shift was selected a priori to approximate the expected timing of the BOLD response. Lag‐specific analyses were restricted to the same subject, run, sampled electrode‐pair connection, frequency‐band, and window keys to ensure that differences across temporal shifts reflected temporal alignment rather than changes in the analytic sample.

Third, to assess sensitivity to dynamic FC window parameters, we repeated the dynamic analyses using shorter 90 s windows with 45 s overlap, longer 150 s windows with 75 s overlap, and non‐overlapping 120 s windows. For each window setting, both icEEG FC and fMRI FC were recomputed using the corresponding window length and overlap. The primary 120 s window with 60 s overlap contained 80 fMRI volumes at TR = 1.5 s and was selected to balance temporal resolution with the reliability of correlation‐based FC estimates. Because a 0.01‐Hz BOLD fluctuation corresponds to a period of approximately 100 s, the primary 120 s window was chosen to provide sufficient temporal coverage for low‐frequency BOLD fluctuations while retaining multiple windows per run. The non‐overlapping 120 s analysis was included to assess whether the dynamic coupling results were influenced by dependence introduced by the 50% overlap between adjacent windows, because overlapping windows share data and may induce serial autocorrelation in windowed FC estimates.

Fourth, to assess whether dynamic cross‐modal coupling varied by anatomical location, the window‐level LMM framework was applied separately to sampled electrode‐pair connections grouped by anatomical region pairs. Bipolar icEEG channels were assigned to retained anatomical bins based on electrode location: LF, left frontal; LMT, left mesial temporal; LI, left insula; RF, right frontal; RMT, right mesial temporal; and RI, right insula. These regions were selected because they were sufficiently represented across the heterogeneous clinical implantation patterns. Channels outside these regions were excluded from this stratified analysis. Because this analysis involved multiple region‐ and frequency‐specific comparisons and was limited by uneven clinically determined electrode coverage, it was treated as exploratory and interpreted descriptively rather than as a confirmatory regional analysis.

Fifth, to test whether icEEG–fMRI coupling depended on anatomical correspondence between modalities, we performed a mismatched‐connection negative‐control analysis. Within each subject and frequency band, fMRI FC values were reassigned to non‐corresponding icEEG electrode‐pair connections, thereby disrupting the anatomical correspondence between electrophysiological and hemodynamic connectivity while preserving the subject‐level and frequency‐specific connectivity distributions. For the dynamic analysis, reassignment was performed within subject, run, window, and frequency band. The same LMM framework was then refit using the mismatched connection pairings. This analysis tested whether the observed coupling was specific to spatially corresponding icEEG–fMRI connections.

Finally, to assess whether the observed icEEG–fMRI coupling was influenced by head motion, we conducted motion‐control robustness analyses. As a first‐level precaution, time windows containing excessive motion, defined as framewise displacement or relative displacement greater than 1.5 mm, were excluded from all analyses. To directly address potential residual motion confounding in static fMRI FC, we repeated the full‐recording static and window‐averaged static FC models with subject‐level mean motion included as an additional covariate. Subject‐level mean motion was computed from run‐level FEAT relative displacement estimates across retained fMRI runs. In addition, to evaluate motion effects in the window‐level analyses, runs were stratified based on run‐level motion severity. For each run, the maximum relative displacement was extracted from the motion‐correction time series generated during fMRI preprocessing in FSL FEAT. Runs were classified as low motion if maximum relative displacement was less than 0.5 mm, and moderate motion if maximum relative displacement was between 0.5 mm and the exclusion threshold. The primary window‐level LMMs were then refit within motion‐defined subgroups to assess whether icEEG–fMRI coupling persisted in low‐motion data.

3.1. Participant Inclusion and Final Analytic Sample

Of the 70 recruited participants, 22 were excluded before the final analysis. Exclusions were due to inability to satisfactorily remove gradient‐switching artifact from the icEEG recordings (n = 5), less than 5 min of continuous fMRI data remaining after motion correction (n = 3), or implantation with strip and/or grid electrodes, either alone or in combination with depth electrodes, rather than depth‐electrode coverage suitable for the present electrode‐pair analysis (n = 14). The final analytic sample therefore included 48 adults with focal epilepsy.

Depth electrodes (Spencer probes SD06R‐SP05X‐000/SD08R‐SP05X‐000/SD10R‐SP05X‐000; Ad‐Tech Medical Instrument Corp. Oak Creek, WI, USA) with 5‐mm inter‐contact spacing were implanted under general anesthesia based on individualized clinical assessment. Only electrodes located in non‐pathological regions, defined as regions without ictal or interictal epileptiform activity, were included in analyses. As expected for clinically determined electrode implantation, regional coverage varied across participants.

Subjects' ages ranged in age from 27 to 72 years (mean ± SD: 44.4 ± 11.3 years) and 22 of the 48 participants were females. On average, 16.6 ± 5.2 non‐epileptogenic regions of interest (ROIs) were analyzed per subject. These ROIs corresponded to electrode contacts in which no IEDs were detected during the recording. This resulted in 143 ± 84 pairwise connections per participant included in the statistical analysis. Detailed demographic and clinical information, together with the effective icEEG‐fMRI analysis duration retained for each participant after motion censoring and exclusion of segments shorter than 5 min, are reported in Table 1.

3.2. Static icEEG Connectivity Is Associated With Spatial Patterns of fMRI FC Across Sampled Electrode‐Pair Connections

Linear mixed‐effects models (LMMs) revealed that mean icEEG FC was positively associated with mean fMRI connectivity across all frequency bands (Figure 4). Connections with stronger electrophysiological coupling consistently showed stronger hemodynamic coupling. When connectivity was estimated over the full retained recording duration, associations were significant in all bands (all $p_{\mathit{FDR}}$ < 0.001; Figure 4A). Effect sizes were largest in the beta and low‐gamma ranges, with comparatively smaller but still significant effects in the delta, theta, alpha, and high‐gamma bands. Full fixed‐effect estimates and 95% confidence intervals for all static models are provided in Table S1.

A complementary analysis using window‐averaged connectivity yielded a similar frequency profile, with systematically larger effect sizes (Figure 4B and Table S1; all $p_{\mathit{FDR}}$ < 0.001). The close agreement between full‐recording static and window‐averaged estimates indicates that time‐averaged windowed connectivity provides a reliable approximation of static connectivity in this dataset, supporting its use as the mean sliding‐window component in subsequent dynamic analyses. Importantly, the positive association between mean icEEG and fMRI FC persisted after adjusting for temporal variability in icEEG connectivity (standard deviation of windowed FC), inter‐electrode distance, age, and sex.

To visualize the spatial organization underlying these statistical associations, representative static connectivity matrices from a single subject are shown in Figure 5. An anatomical rendering of the non‐epileptogenic bipolar contacts and representative pairwise connections is shown alongside the corresponding static fMRI and band‐limited icEEG FC matrices. Bipolar electrodes were ordered by hemisphere, anatomical structure, and orientation to group spatially adjacent contacts. Similar block‐structured patterns were observed across fMRI and icEEG connectivity matrices, particularly in low‐ and intermediate‐frequency bands, in line with the cohort‐level association between electrophysiological and hemodynamic connectivity. Connections between anatomically proximal regions tended to show stronger positive connectivity in both modalities; however, the association between static icEEG FC and fMRI FC remained statistically significant after accounting for inter‐electrode distance. The anatomical panel is provided to visualize the sampled electrode contacts and representative electrode‐pair connections. Statistical inference was based on the full set of eligible sampled electrode‐pair connectivity estimates.

3.3. Window‐Level FC Variance Is Dominated by Within‐Subject Variability

Intraclass correlation (ICC) analysis showed that between‐subject differences accounted for a minority of the total variance in window‐level FC across both fMRI and band‐limited icEEG connectivity estimates (Table 2). For fMRI, the ICC was 0.19 ± 0.05, indicating that 19% of the total variance in windowed fMRI connectivity was attributable to between‐subject differences. Across icEEG frequency bands, ICC values ranged from 0.21 in the theta band to 0.30 in the low‐ and high‐gamma bands. Although ICC increased modestly with frequency in icEEG, values remained well below 0.5 in all bands. Thus, most variance in window‐level connectivity was attributable to within‐subject variation across sampled electrode‐pair connections, runs, windows, and residual variability rather than stable between‐subject differences.

These findings indicate that window‐level connectivity primarily captures within‐subject variation rather than stable trait‐like differences between participants. This variance structure supports the use of window‐level linear mixed‐effects models to examine cross‐modal coupling while accounting for repeated windows, electrode‐pair connections, runs, and subject‐level variability.

3.4. Dynamic fMRI Connectivity Shows Band‐Wise Associations With Window‐Level icEEG Fluctuations

The primary dynamic window‐level LMMs included 216,393 window‐level observations from 4850 sampled electrode‐pair connections across 48 subjects. An omnibus dynamic LMM including a band × window‐specific icEEG FC deviation interaction did not reach conventional significance when compared with a reduced model without the interaction term, although the result suggested possible band‐wise variation in the dynamic icEEG–fMRI coupling effect (likelihood‐ratio test: LR = 9.64, Δdf = 5, $p$ = 0.086). In the primary band‐wise models, window‐specific icEEG FC deviations $\left(\text{ΔicEEG}\right)$ showed significant positive associations with temporally aligned fMRI connectivity in the alpha, beta, and low‐gamma bands. After FDR correction, significant positive dynamic coupling was observed in the alpha band (β = 0.085, 95% CI [0.026, 0.144], $p_{\mathit{FDR}}$ = 0.015) and beta band (β = 0.126, 95% CI [0.038, 0.215], $p_{\mathit{FDR}}$ = 0.015), and low‐gamma band (β = 0.113, 95% CI [0.023, 0.203], $p_{\mathit{FDR}}$ = 0.028). Delta, theta, and high‐gamma bands did not survive FDR correction in the primary band‐wise models (Figure 6A). Together, these results indicate positive window‐level icEEG–fMRI coupling, with the strongest primary band‐wise evidence in alpha, beta, and low‐gamma bands, but without statistically significant omnibus evidence for differences across frequency bands.

In the same window‐level models, the mean sliding‐window icEEG FC term remained positively associated with dynamic fMRI connectivity across all frequency bands (all $p_{\mathit{FDR}}$ < 0.001). This indicates that electrode‐pair connections with stronger average electrophysiological coupling also showed stronger hemodynamic coupling, even after accounting for window‐specific icEEG FC deviations, icEEG FC variability, inter‐electrode distance, age, and sex. The icEEG FC variability term was also associated with dynamic fMRI connectivity across frequency bands. This association was negative in the delta band and positive in theta, alpha, beta, low‐gamma, and high‐gamma bands, suggesting that the overall magnitude of electrophysiological connectivity fluctuations contributed to window‐level coupling patterns in a band‐specific sign pattern.

Allowing subject‐specific slopes for the window‐specific icEEG FC deviation term (ΔicEEG) revealed that while there was a consistent positive group‐level effect, the strength of dynamic icEEG–fMRI coupling varied across individuals. Subject‐specific slope estimates (Figure 6B) showed that although the primary band‐wise group‐level effects were significant in the alpha, beta, and low‐gamma bands, the magnitude and spectral prominence of coupling differed across participants. Some participants also showed negative slopes in specific frequency bands, indicating that, in those individuals, window‐specific icEEG FC deviations and temporally aligned fMRI FC fluctuations tended to vary in opposite directions. Representative examples (Figure 7) illustrate this variability, with one subject showing a greater magnitude of coupling in the beta band and another in the alpha band. These findings indicate that while window‐level cross‐modal dynamic coupling was evident at the group level, its expression varied across individuals and frequency bands.

3.5. Exploratory Spatially Stratified Dynamic icEEG–fMRI Coupling

Exploratory spatially stratified dynamic icEEG–fMRI coupling. To assess spatial variation in cross‐modal coupling, the window‐level LMM analysis was repeated separately for sampled electrode‐pair connections grouped by anatomical region pairs. The resulting 6 × 6 matrices summarize the estimated coupling coefficient (β for ΔicEEG) for each retained region combination, where diagonal entries represent within‐region interactions and off‐diagonal entries represent between‐region interactions. Because this analysis involved multiple region‐ and frequency‐specific comparisons with uneven, clinically determined electrode coverage, regional effects are reported descriptively and were not used to support confirmatory anatomical claims.

Of the 36 possible region–region combinations displayed in the final retained anatomical bins, 28 met predefined inclusion criteria (≥ 6 subjects, ≥ 5 unique electrode pairs, ≥ 500 window‐level observations) and were retained for analysis, whereas 8 combinations (22%) were excluded a priori because of insufficient data coverage. The left and right lateral temporal bins were not retained in the final visualization because they did not meet analysis criteria in any frequency band. This partial coverage reflects clinically driven implantation patterns, with mesial temporal regions more densely sampled than several other regions. The number of contributing subjects, sampled electrode‐pair connections, and window‐level observations for each region pair are reported in Table S6.

The exploratory spatially stratified analysis revealed heterogeneous dynamic icEEG–fMRI coupling patterns across anatomical region pairs and frequency bands (Figure 8). All analyzed cells are colored according to the magnitude and direction of β, with asterisks indicating nominal uncorrected significance levels and n.s. indicating analyzed region pairs that did not pass the nominal threshold. Across all bands and analyzed region pairs, 11 region‐frequency cells showed nominal effects at p < 0.05, whereas 157 analyzed cells did not pass the nominal threshold. Positive estimates indicate that windows with higher‐than‐average icEEG connectivity for a given sampled electrode‐pair connection tended to show stronger temporally aligned fMRI connectivity, whereas negative estimates indicate the opposite pattern.

Across bands, nominal effects were observed in several frontal, mesial temporal, and insular region combinations, with both positive and negative estimates. These patterns should be interpreted cautiously because the analysis was exploratory, unevenly sampled across anatomical regions, and not corrected for all region‐ and frequency‐specific comparisons. Overall, this exploratory analysis suggests that dynamic icEEG–fMRI coupling may vary across implanted anatomical regions and frequency bands, but it should not be interpreted as confirmatory evidence for specific region‐pair effects.

3.6. Sensitivity and Robustness Analyses

Static and window‐averaged icEEG–fMRI coupling remained positive and FDR‐significant across frequency bands when the analysis was repeated using alternative 4‐mm and 8‐mm FWHM smoothing kernels, indicating that average cross‐modal coupling was robust to the choice of spatial smoothing kernel (Table S2A). For the window‐level dynamic models, the 4‐ and 8‐mm analyses were restricted to the same subject, run, electrode‐pair connection, frequency‐band, and window keys used in the primary 6‐mm analysis. The window‐specific icEEG FC deviation term (ΔicEEG) remained positive across all frequency bands in both sensitivity analyses, but the alpha, beta, and low‐gamma effects that survived FDR correction in the primary 6‐mm analysis did not survive correction in the 4‐mm or 8‐mm analyses. Thus, the primary dynamic ΔicEEG effects were more smoothing‐sensitive than the static and average coupling effects (Table S2B).

Temporal‐alignment sensitivity analysis showed that the primary dynamic icEEG–fMRI coupling effects were consistent in direction across fMRI temporal shifts of 0, 3, 6, and 9 s (Table S3). After restricting all lag‐specific analyses to the same subject, run, sampled electrode‐pair connection, frequency‐band, and window keys, ΔicEEG effects were positive across all frequency bands and all temporal shifts. The strongest and most consistent statistical evidence was observed in the alpha, beta, and low‐gamma bands, which remained FDR‐significant at the 3, 6, and 9 s shifts, whereas the corresponding 0 s effects were positive but did not survive FDR correction. Delta, theta, and high‐gamma effects were also positive across temporal shifts but did not survive FDR correction. These results indicate that the direction of the dynamic icEEG–fMRI coupling effect was not dependent on a single assumed temporal alignment, although statistical significance was strongest when fMRI windows were shifted to account for the expected BOLD delay.

Window‐parameter sensitivity analysis showed that the direction of the primary dynamic icEEG–fMRI coupling effects was consistent across alternative window definitions (Table S4). ΔicEEG effects were positive across all frequency bands when dynamic FC was recomputed using 90 s windows with 45 s overlap, the primary 120 s windows with 60 s overlap, 150 s windows with 75 s overlap, and non‐overlapping 120 s windows. The strongest and most consistent statistical evidence was observed in the alpha, beta, and low‐gamma bands, which survived FDR correction in the 90 s, primary 120 s, and 150 s overlapping‐window analyses. In the non‐overlapping 120 s analysis, effects remained positive across all bands but did not survive FDR correction. Theta also reached FDR significance only in the 150 s window analysis, whereas delta and high‐gamma effects did not survive FDR correction under any window setting. Overall, these results indicate that the direction of dynamic icEEG–fMRI coupling was not dependent on a single window length or overlap choice, although statistical significance was strongest for overlapping‐window analyses.

In the mismatched electrode‐pair negative‐control analysis, disrupting the anatomical correspondence between icEEG and fMRI sampled electrode‐pair connections eliminated the static coupling effects and strongly attenuated the dynamic effects (Table S5). Static mismatched‐connection associations were near zero and did not survive FDR correction in any frequency band. Dynamic mismatched‐connection effects were also substantially smaller than the primary anatomically matched effects and did not survive FDR correction in any frequency band. Together, these findings indicate that the observed static and dynamic icEEG–fMRI coupling depended on anatomical correspondence between sampled icEEG and fMRI electrode‐pair connections.

To evaluate whether dynamic icEEG–fMRI coupling was influenced by head motion, we repeated the primary window‐level analyses separately within low‐motion and moderate‐motion runs (Figure 9). Across frequency bands, coupling coefficients remained positive and comparable in magnitude within the low‐motion subgroup. Confidence intervals substantially overlapped between motion strata. Importantly, coupling estimates were numerically stronger under low‐motion conditions, with no evidence of attenuation relative to higher‐motion segments. These findings indicate that the observed cross‐modal coupling effects are unlikely to be driven by residual head motion artifacts and may in fact be more robust when motion is minimized.

In an additional motion‐adjusted sensitivity analysis, the static icEEG–fMRI coupling effects remained positive and FDR‐significant across all frequency bands after including subject‐level mean motion as an additional covariate for both full‐recording static FC and window‐averaged static FC models (Table S7). These results indicate that the static coupling findings were not explained by residual head motion.

4.1. Static Electrophysiological Coupling Is Associated With Spatial Variation in BOLD Connectivity

Across all frequency bands, time‐averaged icEEG connectivity was significantly associated with static fMRI connectivity. Connections with stronger electrophysiological coupling exhibited stronger haemodynamic coupling, even after controlling for inter‐electrode distance, age, sex, and residual head motion in sensitivity analyses. This broadband correspondence indicates that spatial variation in resting‐state BOLD connectivity is closely related to spatial variation in electrophysiological coupling across sampled electrode‐pair connections. Importantly, the association was preserved whether connectivity was computed across the full retained recording or obtained by averaging windowed estimates, indicating that the spatial correspondence was not simply a consequence of the windowing procedure. Rather, haemodynamic connectivity appears to reflect electrophysiological coupling accumulated over time across ongoing activity (Deligianni et al. 2014; Magri et al. 2012; He et al. 2008; De Pasquale et al. 2010; Pfurtscheller and Lopes Da Silva 1999).

4.2. Window‐Level Electrophysiological Fluctuations Show Band‐Wise Associations With Dynamic BOLD Connectivity

Beyond static coupling, window‐specific deviations in icEEG connectivity from the mean sliding‐window icEEG FC were associated with temporally aligned fMRI connectivity in the primary band‐wise models. In these models, significant positive dynamic coupling was observed in the alpha, beta, and low‐gamma bands after FDR correction, whereas delta, theta, and high‐gamma bands did not survive correction. However, an omnibus band × ΔicEEG interaction did not reach conventional significance, although it suggested possible band‐wise variation in the dynamic coupling effect. Therefore, the band‐wise findings should not be interpreted as statistically confirmed evidence for differences between frequency bands. Rather, the results indicate that dynamic BOLD FC showed the strongest primary band‐wise evidence of association with window‐specific icEEG FC deviations in the alpha, beta, and low‐gamma bands. These effects were estimated after accounting for mean sliding‐window icEEG FC, overall icEEG FC variability, inter‐electrode distance, age, sex, and subject‐specific variability in the ΔicEEG effect.

Interestingly, the presence of negative subject‐specific slopes in a few subjects indicates that dynamic icEEG–fMRI coupling was not uniformly positive across individuals. In these participants and frequency bands, window‐specific icEEG FC deviations and temporally aligned fMRI FC tended to vary in opposite directions. These inverse associations may reflect heterogeneity in anatomical sampling, inter‐individual variability in neurovascular coupling, differences in vigilance or brain state, medication and disease‐related effects, or subject‐specific signal‐to‐noise differences. Because the present study was not designed to classify participants based on slope direction, these negative slopes should be interpreted cautiously. Future work with larger samples and more systematic anatomical sampling will be needed to determine whether inverse dynamic coupling reflects meaningful physiological subtypes or measurement variability.

Prior invasive electrophysiological studies in non‐human primates established that local field potential activity is closely related to BOLD responses, providing foundational evidence for the neuronal basis of the BOLD signal (Logothetis et al. 2001). However, those studies primarily addressed local electrophysiological activity in relation to regional BOLD amplitude, whereas the present study examined inter‐regional connectivity in simultaneous icEEG and fMRI. In addition, our high‐gamma analysis was restricted to 51–100 Hz, which should be interpreted as a restricted high‐gamma range rather than broadband high‐gamma activity extending to 180–200 Hz or higher. We therefore interpret the absence of robust dynamic coupling in the 51–100 Hz window‐specific icEEG FC deviation term cautiously. This finding does not rule out a relationship between broadband high‐gamma activity and BOLD responses but indicates that, within the frequency range and sampled electrode‐pair connectivity framework analyzed here, dynamic icEEG–fMRI coupling showed the strongest primary band‐wise evidence for dynamic icEEG–fMRI coupling observed in the alpha, beta, and low‐gamma bands.

The alpha‐band result should be distinguished from prior reports of alpha–BOLD anticorrelation (Goldman et al. 2002; Laufs et al. 2003). These simultaneous scalp EEG–fMRI studies were conducted during eyes‐closed rest and primarily examined regional alpha power in relation to local BOLD amplitude. Although both studies reported negative associations between alpha power and cortical BOLD signal, positive alpha‐related BOLD changes have also been reported in thalamic and insular regions (Goldman et al. 2002). In contrast, the present study examined inter‐regional icEEG connectivity in relation to inter‐regional fMRI connectivity across sampled electrode‐pair connections. Therefore, the positive alpha‐band icEEG–fMRI connectivity association observed here does not contradict prior alpha–BOLD anticorrelation findings, but instead suggests that alpha–BOLD relationships may depend on the electrophysiological feature measured, the spatial scale of analysis, and whether local amplitude or inter‐regional connectivity is examined. Because cognitive state, vigilance, arousal, and behavioral task demands were not directly measured, we do not infer specific functional roles for the alpha‐band coupling observed here.

Delta‐ and theta‐band window‐specific icEEG FC deviations did not show robust dynamic coupling with temporally aligned fMRI FC in the primary band‐wise models. We do not interpret this as evidence that low‐frequency activity is unrelated to BOLD dynamics. Rather, sensitivity to low‐frequency coupling may depend on features of the present windowed FC framework, including the use of correlation‐based connectivity estimates over finite time windows, the focus on electrode‐pair band‐limited power‐correlation FC rather than local power–BOLD amplitude coupling, and the use of fixed temporal alignments between icEEG and fMRI. Low‐frequency effects may also occur over slower timescales than those captured by the present window‐level deviation analysis. Thus, the non‐significant delta/theta findings should be interpreted as specific to the present windowed FC framework rather than as evidence against low‐frequency neural contributions to BOLD activity.

The mismatched electrode‐pair negative‐control analysis further supports the anatomical specificity of the observed cross‐modal associations. When the correspondence between sampled icEEG and fMRI electrode‐pair connections was disrupted within subjects, static coupling effects were eliminated and dynamic ΔicEEG effects were strongly attenuated. This indicates that the observed icEEG–fMRI coupling was not simply driven by subject‐level differences or overall connectivity distributions, but instead depended on anatomically corresponding connectivity estimates across modalities.

4.3. Exploratory Spatial Variation of Dynamic Coupling Across Implanted Regions

The spatially stratified analysis suggested that dynamic icEEG–fMRI coupling may vary across the retained clinically sampled anatomical region pairs, including frontal, mesial temporal, and insular combinations. However, because electrode coverage was clinically determined, uneven across participants, and the nominal regional analyses were not corrected across all region‐ and frequency‐specific comparisons, these findings should be interpreted descriptively rather than as confirmatory evidence for specific anatomical effects. Positive and negative regional estimates indicate that the direction of coupling between window‐specific icEEG FC deviations and temporally aligned fMRI connectivity may differ across sampled region pairs, but the biological significance of these effects remains uncertain. Future studies with systematic anatomical sampling and appropriate correction for regional multiple comparisons will be needed to determine whether these spatial patterns reflect reproducible biological effects.

4.4. Motion Robustness

Dynamic coupling estimates were preserved when analyses were restricted to low‐motion runs. Overlapping confidence intervals across motion strata and the absence of systematic attenuation suggest that the observed cross‐modal associations were unlikely to be driven by residual head motion artifacts. In addition, static icEEG–fMRI coupling remained positive and FDR‐significant across all frequency bands after including subject‐level mean motion as an additional covariate, further supporting the robustness of the static coupling findings.

4.5. Relation to Prior Simultaneous icEEG–fMRI Studies

Simultaneous icEEG–fMRI datasets remain rare, limiting direct comparison. Ridley et al. (2017) reported cross‐modal correlations between icEEG and BOLD connectivity, strongest in lower‐frequency bands and non‐epileptogenic cortex. Our results are broadly consistent with this finding, demonstrating robust static associations between electrophysiological and haemodynamic connectivity across frequency bands. The present study extends this prior work by separating average sliding‐window electrophysiological connectivity from window‐level deviations around that average, allowing static and dynamic cross‐modal associations to be examined within the same modeling framework.

Mostame et al. (2025) proposed that large‐scale connectome dynamics operate as a multiplex system across timescales, showing that spatially similar icEEG‐ and fMRI‐derived states recur asynchronously. Extending this framework, we quantified cross‐modal coupling at the level of individual electrode‐pair connections rather than global state recurrence. Our findings suggest that temporal asynchrony at the global state level does not preclude window‐level coupling at the sampled electrode‐pair connection level.

A key distinction of the present study is scale. Prior simultaneous icEEG–fMRI investigations were based on small cohorts, limiting population‐level inference. Leveraging 48 subjects and hierarchical mixed‐effects modeling, we provide a population‐level characterization of cross‐modal coupling across individuals and clinically sampled anatomical contexts. Restricting analyses to non‐epileptogenic electrodes further reduced the influence of local epileptogenic activity on the observed cross‐modal associations.