A multiplex of connectome trajectories enables several connectivity patterns in parallel

Abstract

Complex brain function comprises a multitude of neural operations in parallel and often at different speeds. Each of these operations is carried out across a network of distributed brain regions. How multiple distributed processes are facilitated in parallel is largely unknown. We postulate that such processing relies on a multiplex of dynamic network patterns emerging in parallel but from different functional connectivity (FC) timescales. Given the dominance of inherently slow fMRI in network science, it is unknown whether the brain leverages such multi-timescale network dynamics. We studied FC dynamics concurrently across a breadth of timescales (from infraslow to γ-range) in rare, simultaneously recorded intracranial EEG and fMRI in humans, and source-localized scalp EEG-fMRI data in humans. We examined spatial and temporal convergence of connectome trajectories across timescales. 'Spatial convergence' refers to spatially similar EEG and fMRI connectome patterns, while 'temporal convergence' signifies the more specific case of spatial convergence at corresponding timepoints in EEG and fMRI. We observed spatial convergence but temporal divergence across FC timescales; connectome states (recurrent FC patterns) with partial spatial similarity were found in fMRI and all EEG frequency bands, but these occurred asynchronously across FC timescales. Our findings suggest that hemodynamic and frequency-specific electrophysiological signals, while involving similar large-scale networks, represent functionally distinct connectome trajectories that operate at different FC speeds and in parallel. This multiplex is poised to enable concurrent connectivity across multiple sets of brain regions independently.

Keywords: concurrent; connectome dynamics; fMRI; functinoal connectome; human; intracranial EEG; multiplex; neuroscience.

Figure 1.. Three potential scenarios for convergence/divergence of connectome trajectories across hemodynamics and electrophysiology.

Dynamics of connectome reconfigurations (i.e. time-varying spatial patterns of large-scale connectivity) are shown as trajectories (functional connectivity, FC pattern sequences) in a presumed 3D state space, starting from the diamond sign and ending at the arrowhead. Hemodynamics- and electrophysiology-derived connectome dynamics in each scenario are shown in blue and red, respectively. Scenario (I): both dynamic trajectories are driven by the same connectivity processes; are spatially and temporally convergent (at all times). Scenario (II): both dynamic trajectories span the same portion of the state space, while each traversing a different timecourse; are spatially convergent but temporally divergent. Scenario (III): each of the dynamic trajectories spans non-overlapping parts of the state space with different timecourses; are spatially and temporally divergent (at all times). Note that scenarios I and III represent border cases extracted from a theoretical convergence/divergence gradient (bottom of the figure). Spatial convergence does not imply a perfect spatial match (at the same or different time points). Rather, spatial convergence –if observed– is likely to be partial, given that the cross-modal similarity of static connectomes is known to be moderate (Betzel et al., 2019; Wirsich et al., 2021). Furthermore, note that electrophysiology-derived connectivity is simplified for illustration purposes and may comprise several trajectories at different frequency bands.

Figure 1—figure supplement 1.. Electrode locations for each intracranial EEG (iEEG)-fMRI subject shown in standard MNI space.

Only the ECoG and depth electrode contacts that were included in the analyses are shown (black dots). In particular, electrodes that were directly or indirectly related to the cause of seizures (spike analysis described in Ridley et al., 2017), contained epileptiform activity, or were embedded in white matter were excluded.

Figure 2.. Cross-modal recurrence plots (CRPs) quantify synchronous and asynchronous spatial convergence across data modalities.

(A) CRP between fMRI-FC and iEEG-FC_Amp shown for different electrophysiological bands (δ to γ) for a single subject. Note that synchronous observations (after accounting for a 6 s hemodynamic delay) fall on the diagonal of the CRP. The CRP in each frequency band is binarized, with correlation values passing the threshold if they exceed a null model generated from spatially phase-permuted FC matrices (Benjamini-Hochberg FDR-corrected for the number of all pairs of timepoints). (B) Multi-frequency CRP constructed for the same subject by overlaying the CRPs across all electrophysiological frequency bands. Pairs of timepoints with significant spatial correlation across fMRI and intracranial EEG (iEEG) functional connectivity (FC) are indicated with a different color for each frequency band. (C) Schematic views of expected CRP patterns corresponding to each of the three scenarios introduced in Figure 1. Significant CRP entries are marked black. The density of black bins depicts the extent of spatial convergence, while their prominence along the diagonal (measured as on-/off-diagonal ratio) translates to the extent of temporal convergence. In the following, the density and on-/off-diagonal ratio of significant entries in the multi-frequency CRP will be used to adjudicate the likelihood of the scenarios depicted in Figure 1 (see section ‘Validation of on-/off-diagonal ratio metric’ in Appendix 1).

Figure 3.. The spatial convergence across fMRI and frequency-specific intracranial EEG (iEEG) connectomes occurs asynchronously.

(A) Effect size of cross-modal correlation: Frequency bands are color-coded using the same set of colors from Figure 2. Within each frequency band, each pair of opaque and transparent violin plots corresponds to one subject. Opaque (or transparent) plots represent the distribution of correlation values in significant (or insignificant) multi-frequency cross-modal recurrence plot (CRP) entries. The horizontal black lines in each violin plot represent the mean value of the distribution. Across all frequency bands, the effect size of spatial correlation between fMRI and iEEG-FC_Amp connectome frames at significant epochs is considerably larger than in insignificant epochs. (B) Total density of significant epochs: Horizontal axis shows the subject number (rightmost column shows group average), while vertical axis shows percentage of significant epochs across all multi-frequency CRP entries. Overall, 20% of multi-frequency CRP entries averaged over subjects passed significance testing, indicating spatial convergence between fMRI and iEEG connectomes and speaking against scenario III. (C) On-/off-diagonal ratio of significant epochs calculated across a range of diagonal shifts in CRP. The horizontal axis shows a range of lead/lags by 5 TRs at which fMRI and iEEG were assumed to be aligned (i.e. ‘on-diagonal’) (see Figure 2). Recall that zero lag in CRP corresponds to original preprocessed data where fMRI is shifted 6 s to account for hemodynamic delay. Gray dots show the ratio for each individual and the horizontal lines represent the group-average ratios at each lag. Overall, the ratios are very close to 1 and do not exceed data-driven null distribution. This observation speaks in favor of scenario II where fMRI and iEEG connectome reconfigurations follow distinct temporal trajectories while manifesting spatially similar connectome patterns. (D) Degree of overlap between pairs of frequency-specific binarized CRPs (as shown in Figure 2A), expressed as group-average of the Jaccard index. The low values (cf. possible index range of 0–1) demonstrate that timepoint pairs of high cross-modal spatial correlation are largely distinct across different electrophysiological frequency bands. In summary, the observed cross-modal spatial convergence at frequency-specific and asynchronous epochs suggests that connectome dynamics across different functional connectivity (FC) timescales traverse spatially similar patterns asynchronously (Scenario II; Figure 5).

Figure 3—figure supplement 1.. Static intracranial EEG (iEEG)/EEG connectome as an attractor for the cross-modal convergence.

(A) Quantitative assessment of the emergence of cross-modal recurrence plot (CRP) stripes in both intracranial EEG-fMRI (top) and scalp EEG-fMRI (bottom) datasets. The skewness of the cumulative number of significant elements (or ‘Peakiness’) along each CRP axis can capture structured patterns, specifically stripes in horizontal and vertical directions. For fMRI-iEEG data, the group-average skewness along fMRI axis (Skewness = 1.60 ± 0.65) was significantly higher than along the multi-frequency iEEG-FC_Amp axis (Skewness = 0.49 ± 0.36; paired t-test between axes and across subjects: t₈=4.05 p=0.0018). A similar difference was observed for the source-localized EEG-fMRI dataset (t₂₅=10.83; p=3.14e-11). Importantly, the higher skewness along the horizontal (fMRI) axis than the vertical (iEEG/EEG) axis was observed in individuals (7 out of 9 and 25 out of 26 subjects, respectively in iEEG-fMRI and EEG-fMRI datasets). Thus, the spatial similarity during the stripes occurs largely because of spatial characteristics that are present in electrophysiology in a stable manner over time, and are shared across frequency bands (since stripes were stronger when all frequency-specific CRPs were combined). Therefore, the time-averaged electrophysiology-derived connectome may serve as an attractor for the cross-modal convergence of connectome dynamics. (B) The schematic view of the approach for testing whether the time-averaged electrophysiology-derived connectome serves as an attractor for the cross-modal convergence. The multi-frequency CRP (left side), the accumulated number of significant entries in the CRP along vertical axis (middle), and Pseudo-CRP (right) are depicted for a representative subject of the source-localized EEG-fMRI data. The Pseudo-CRP was defined as the correlation of fMRI connectomes to the time-averaged EEG connectome. (C) In each subject, temporal correlation between the accumulative counts of significant entries in the multi-frequency CRP and the Pseudo-CRP (blue dots) were compared to a null model derived from randomized Pseudo-CRPs (gray dot clouds; R=1000). We found a strong correlation across all subjects of the intracranial EEG-fMRI (real data: r=0.73 ± 0.09; for null: r=0.01 ± 0.11; p<10^–3 for all subjects) and source-localized EEG-fMRI (real data: r=0.88 ± 0.03; for null: r=0.01 ± 0.12; p<10^–3 for all subjects) datasets. Significant temporal correlation between the CRP stripiness timecourse and the Pseudo-CRP timecourse suggests that the time-averaged EEG connectome serves as an attractor for cross-modal correlations of fMRI and EEG connectome dynamics.

Figure 3—figure supplement 2.. Temporal association between connection-wise dynamic functional connectivity (FC) changes of fMRI and intracranial EEG (iEEG).

(A) Significance rate of the temporal association of dynamic FC changes is extracted using the mutual information index, and then comparing with a phase-permuted null model (50 iterations; FDR-corrected for the number of edges within each subject). The significance rate of the temporal association (vertical axis) is reported individually (first five sets of columns; ‘Delta’ to ‘Gamma’) and cumulatively across all frequency bands (last set of columns; ‘All-frequencies’). Each box indicates first, second, and third quartiles of the data, while whiskers extend to the last data sample that is not considered an outlier (marked with red plus). Samples that are 1.5 times the interquartile range larger (or smaller) than the third (or first) quartile are considered outliers. A considerable proportion of the connections show temporal association between fMRI and iEEG FC changes at the connection level. Importantly, this association cumulatively increases when frequency bands in electrophysiology are pooled together. This observation speaks to partially shared multi-frequency neurophysiological processes driving portions of both fMRI- and electrophysiology-derived FC at the connection level.

Figure 3—figure supplement 3.. Replication of the results in the scalp EEG dataset after correction for volume conduction effects using source-orthogonalized signals for both amplitude- and phase-coupling measures (similar to Figure 4 in the main text).

(A) Overall, ~15% of the multi-frequency cross-modal recurrence plot (CRP) entries showed significant correlations between fMRI-FC and source-localized EEG-FC_Amp (~12% for source-localized EEG-FC_Phase), speaking to a temporally sparse spatial convergence between fMRI and EEG connectome dynamics. (B) The on-/off-diagonal ratio (at zero diagonal shift in CRP) across subjects was very close to 1 and did not exceed null data for both amplitude coupling (real ratio = 0.99 ± 0.15; mean null ratio = 1.00 ± 0.01; t₂₅=–0.11; p=0.54; BF₀₁=5.2) and phase coupling (real ratio = 0.97 ± 0.16; mean null ratio = 1.00 ± 0.02; t₂₆=–0.98; p=0.83; BF₀₁=8.8). These findings were replicated across different diagonal shifts in CRP for both amplitude coupling (t₂₅=−0.94 ± 0.90; FDR-corrected p>0.99; mean ratio = 0.99 ± 0.02, mean null ratio = 1.01 ± 0.01, BF₀₁=8.7 ± 3.4) and phase coupling (t₂₅=−0.05 ± 0.68; FDR-corrected p>0.85; mean ratio = 1.01 ± 0.02, mean null ratio = 1.01 ± 0.01, BF₀₁=5.2 ± 2.4). These observations, again, suggest temporally independent connectivity processes between the two modalities. (C) The cross-frequency overlap of significant bins in the multi-frequency CRP was negligible (Jaccard index <0.03 for fMRI-FC to source-localized EEG-FC_Amp and <0.02 for fMRI-FC to source-localized EEG-FC_Phase), confirming the observed multi-frequency nature of the association between fMRI and EEG connectome dynamics.

Figure 3—figure supplement 4.. Replication of major findings using phase coupling as a measure of electrophysiology connectivity in intracranial EEG and source-localized scalp EEG datasets.

(A) When using phase coupling measure, respectively in the intracranial and scalp EEG data, ~18 and ~30% of the multi-frequency cross-modal recurrence plot (CRP) entries showed significant correlations between fMRI-FC and iEEG-FC_Phase, which speaks against scenario III. (B) The on-/off-diagonal ratio (at zero diagonal shift in CRP) in each subject did not exceed the subject-specific null model (p<0.11). At the group-level, the ratio was very close to 1 and did not exceed surrogate data for both iEEG (real ratio = 1.00 ± 0.17; mean null ratio = 1.00 ± 0.02; t₈=0.10; p=0.46; BF₀₁=2.9) and scalp EEG (real ratio = 0.99 ± 0.08; mean null ratio = 1.00 ± 0.01; t₂₆=–0.19; p=0.57; BF₀₁=5.5) data. These findings were replicated across different diagonal shifts in CRP for both iEEG (t₈=−0.66 ± 1.53; FDR-corrected p>0.99; mean ratio = 1.00 ± 0.08, mean null ratio = 1.02 ± 0.01, BF₀₁=4.4 ± 2.7) and scalp EEG (t₂₅=−1.21 ± 1.26; FDR-corrected p>0.99; mean ratio = 0.99 ± 0.02, mean null ratio = 1.01 ± 0.01, BF₀₁=10.0 ± 4.6) datasets. These results speak to the presence of an asynchronous convergence of connectome configurations across modalities (Scenario II). (C) The cross-frequency overlap of significant bins in the multi-frequency CRP was negligible (Jaccard index <0.12 for fMRI-FC to iEEG-FC_Phase and <0.13 for fMRI-FC to scalp EEG-FC_Phase), confirming the observed multi-frequency nature of the association between fMRI and EEG connectome dynamics.

Figure 3—figure supplement 5.. The on-/off-diagonal ratio is highly replicable when using a window size of 500 ms compared to original findings with a window size of one TR.

On-/off-diagonal ratio is shown across lags for both analyses: 3 s window length (blue) and 500 ms window length (red). Each bar shows the mean across subjects, while the whiskers show the corresponding standard deviations.

Figure 3—figure supplement 6.. cross-modal recurrence plot (CRP) calculated for simulated data with two different cases.

Left: Two identical state time courses resulting in high diagonal prominence, in line with Scenario I, where both spatial and temporal convergence are present. Right: A state time course and its randomized version, showing dispersed co-occurrences and reduced diagonal prominence, representing a case with spatial convergence but no temporal convergence (Scenario II).

Figure 3—figure supplement 7.. On-/off-diagonal ratio across varying quantization levels (N) for simulated data with randomized signals (Scenario II) and identical signals (Scenario I).

The orange bars show the ratio for the simulated data with randomized signals, which stay close to 1 across all Ν values. This observation confirms the absence of temporal synchronization in line with Scenario II. The blue bars show the ratio for the simulated data with identical signals, increasing consistently with N. This observation indicates strong temporal synchrony in line with Scenario I. These results validate the metric’s ability to distinguish between Scenario I and Scenario II.

Figure 4.. Findings generalize to the whole-brain connectome in concurrently recorded source-localized scalp EEG-fMRI.

(A) Total significance rate of the cross-modal spatial similarity passing significance testing in the multi-frequency cross-modal recurrence plot (CRP), shown across all subjects (N=26; analogous to Figure 3B). (B) The on-/off-diagonal ratio of the rate of spatially correlated epochs in the multi-frequency CRP for different diagonal shifts (analogous to Figure 3C). Note that lag of zero corresponds to the original preprocessed data where fMRI has been shifted 6 s to account for the hemodynamic delay. (C) Overlap of spatially correlated epochs across pairs of electrophysiological frequencies in the multi-frequency CRP (averaged over subjects; analogous to Figure 3D). Equivalent findings are presented for EEG phase coupling connectomes as a supplemental analysis (see Appendix 1 and Figure 3—figure supplement 4).

Figure 5.. Emergence of canonical intrinsic connectivity networks (ICNs) across the breadth of functional connectivity timescales.

In a sample subject, overall connectivity strength within the Visual network (top) and the Default mode network (bottom) are separately shown for all functional connectivity (FC) timescales during the whole recording period (fMRI-FC and EEG-FC_Amp, averaged across all region pairs within the respective network). Left: canonical Visual and Default mode networks (Yeo et al., 2011) rendered on a cortical surface. Center: The timeseries show the continuous changes of FC strength (z-scored for better visualization) for all timescales sorted from slow to fast in ascending order, using the same color coding as in Figure 2. Thresholded strength of each network (z>2) is shown for all timescale in a similar fashion on the bottom section of the plot to emphasize periods of particularly high connectivity. The gray vertical bar marks a sample timepoint where the two networks emerge concurrently at two different timescales: the Visual network in the EEG-theta band and the Default mode network in fMRI. This observation demonstrates the co-occurrence of multiple distinct networks at different FC timescales at one given timepoint. Right: The matrices show the whole-brain spatial FC patterns averaged over the thresholded timepoints (z>2) of the respective timescale. The cross-timescale spatial similarity is visually apparent among the top three (respectively bottom three) matrices. The complete set of FC timescales is shown in Figure 5—figure supplement 1. In line with scenario II, these spatially similar FC patterns (dominated by high FC within the Visual or Default mode network, respectively) occur asynchronously across timescales. Please note that this figure aims at making the cross-timescale spatial similarity of connectomes visually accessible and does not reflect a quantitative analysis. For the statistical approach and analyses of this study, see Figures 2—4.

Figure 5—figure supplement 1.. The matrices are the same as shown on the right side of Figure 5 in the main text, but for all timescales (fMRI and δ-γ in EEG).

Matrices represent the average spatial pattern manifested in the connectome dynamics of a particular timescale (fMRI and δ-γ in EEG), only including frames that correlate to the Visual network (top) or Default mode network (bottom) binary mask above an arbitrary threshold. Although we acknowledge the circularity of this approach, note that the purpose of this figure is to make the cross-timescale spatial similarity of connectomes visually accessible and does not reflect a quantitative analysis. For the statistical approach and analyses of this study see Figures 2—4 in the main text.

Figure 6.. fMRI and EEG connectome convergence arises from distinct recurrent states.

(A) Left: Group-level static connectomes of fMRI (top) and EEG (bottom) are illustrated. Right: Five dissociable connectome states at the group level were identified for fMRI and EEG (depicted for β band as an example). Red and blue colors represent positive and negative connectivity values, respectively. (B) Similar to A but for a sample surrogate dataset with randomly phase-permuted connectivity timecourses. Note that states in the surrogate data are highly similar to one another and the static connectome. (C) The vertical axis shows the extent to which the detected connectome states within a functional connectivity (FC) timescale are dissociable. The state dissociation was calculated by the similarity between the proximity matrix of the empirical data and the proximity matrix extracted from cluster labels (a binary proximity matrix representing complete dissociation). The horizontal axis shows timescales (fMRI and EEG frequency bands) on which the clustering was performed. Filled circles show the real effect while the transparent circles show surrogate data. In every timescale, K-means states in the real data were statistically more dissociable than every single surrogate sample (z=103.6, 128.0, 91.2, 35.5, 5.5, 3.7 for fMRI and δ- to γ-bands of EEG, respectively). These observations speak to the presence of distinct recurrent patterns among the converging connectome patterns across timescales.

Figure 6—figure supplement 1.. Replication of clustering analysis with number of states equal to seven.

(A and B) show the static connectomes and the seven states in both real and null data, analogous to Figure 6A and B in the main text C. The vertical axis shows the extent to which the detected connectome states within a timescale are dissociable (measured by proximity matrix). The horizontal axis shows timescales (fMRI and EEG frequency bands) on which the clustering was performed. Filled circles show the real effect while the transparent circles show surrogate data. Connectome states with K=7 were considerably more distinct from one another compared to every single surrogate data (z = 95.1, 94.9, 54.1, 22.1, 5.1, 3.0 for fMRI, and δ to γ bands of EEG).

Figure 7.. Multiplex of dynamic connectome trajectories across timescales of functional connectivity captured by fMRI and electrophysiology.

Our findings indicate that connectome dynamics across fMRI and each frequency band of electrophysiology are spatially convergent (comprise sets of similar spatial states) but temporally divergent (the states occur at different times). Co-existence of such temporally divergent connectome trajectories across functional connectivity (FC) timescales is consistent with a multiplex of connectome patterns at any particular time point. (A) Each color-coded trajectory represents the dynamic spatial reconfigurations of connectome patterns of a particular timescale in state space, within a hypothetical time period. The dots marked on each trajectory represent the same instant in time (e.g. t₁₀), highlighting the different position that each trajectory occupies in state space at this instant. Note that trajectories are partially occupying the same area of the state space (c.f. scenario II in Figure 1). (B) An example of the co-existence of connectome patterns across different timescales embedded within a multiplex network at the exemplary timepoint t₁₀. Functionally, this multiplex of connectome configurations may serve as parallel ‘communication channels’ that maximize information flow between brain regions, akin to the multitude of independent information streams multiplexed into a single radio or television transmission medium.

Figure 7—figure supplement 1.. Alternative sources of contributions to the cross-modal similarity timecourse.

(A) Assessing the temporal correlation between the cross-modal similarity timecourse and the frame-by-frame static connectome ‘prominence’ timecourse of the connectome configurations. Left panel shows results for the intracranial EEG-fMRI dataset while the right panel shows results for the scalp EEG-fMRI dataset. For each frequency band (column pairs), the left column represents comparison to the static fMRI connectome, and the right column the comparison to the intracranial EEG (iEEG)/EEG static connectome. Colored dots show real temporal correlations and gray dots show corresponding null data pooled across subjects. Same color code as in Figure 2 is used for the frequency bands. (B) Same as A but for assessing the temporal correlation of the cross-modal similarity timecourse and the overall (root mean square of all connections) FC strength. (C) Assessing temporal overlap of the cross-modal similarity timecourse and head motion across both datasets, five frequency bands, and all subjects (analogous to A&B). (D) Same as C, except for assessing the temporal overlap between the binarized cross-modal similarity timecourse and epileptiform activity in the iEEG-fMRI dataset. Altogether, these findings suggest that our main findings are not associated with any of the mentioned artefactual or confounding factors.