The coupling between the spatial and temporal scales of neural processes revealed by a joint time-vertex connectome spectral analysis

the relationship be-tween spatially localized or distributed neural processes on one side and their respective temporal frequency bands in source-reconstructed M/EEG signals. We explore this approach on two different datasets, an auditory steady-state response (ASSR) and a visual grating task. Our results suggest that different information processing mechanisms are carried out at different frequency bands: while spatially distributed activity (which may also be interpreted as integration) specifically occurs at low temporal frequencies (alpha and theta) and low graph spatial frequencies, localized electrical activity (i


Introduction
Brain oscillations are produced by the coordinated activity of large groups of neurons and are interpreted as intermittent changes in neuronal excitability (Buzsáki and Draguhn, 2004). Human brain oscillations can be measured at the scalp by electroencephalography (EEG). Traditional approaches to EEG analysis involve the decomposition of the EEG signal into functionally distinct temporal scales or frequency bands, each of which has different temporal and topographical characteristics. In adults, typical frequency bands and their approximate spectral boundaries are delta (1-3 Hz), theta (4-7 Hz), alpha (8-12 Hz), beta (13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30), and gamma (30-100 Hz) (Buzsáki and Draguhn, 2004;Lakatos et al., 2005). Distinct brain rhythms are thought to support different modes of information processing. Theta oscillations have been proposed to serve the attentional sampling of the incoming sensory stimuli (Fries, 2015). Alpha and beta oscillations are thought to mediate top-down attentional mechanisms (Fries, 2015), while gamma-band oscillations are thought to reflect the propagation of incoming stimuli in the respective primary sensory cortices (Fries, 2015). In addition, neurophysiological studies reveal the specificity of cortical layers in the visual system. Theta and gamma rhythms originate from the supragranular layers and serve as feedforward signaling, while feedback influences are carried by alpha and beta band oscillations which arise from infragranular layers (Bastos et al., 2015;van Kerkoerle et al., 2014).
For these reasons, brain oscillations are thought to be an important feature of large-scale network dynamics, with the potential to link specific patterns of neural activity to behavior (Buzsáki and Draguhn, 2004). Recently, it has been shown that the oscillatory activity in distinct frequency bands is differentially constrained by the anatomical architecture of the white matter -the connectome -in primates (Vezoli et al., 2021). However, how distinct brain rhythms support different modes of information processing at the different scales in the human brain remains unknown.
The term integration has to be understood as a global communication mechanism that combines and distributes neural information throughout the brain, ensuring large-scale coherent functional activity. Segregation entails the notion of functional specialization, in which information is relayed predominantly inside localized brain regions involved in a specialized function. Integration and segregation have been described as the two fundamental mechanisms of information processing to enable proper sensory processing (Shine, 2019), coordinated behavior and flexible cognition (Deco et al., 2015). Several neuroimaging methods allow the study of segregation and integration (Deco et al., 2015): at the microscopic level, neural assemblies are defined as distributed networks of neurons, that are linked by reciprocal connections that can be direct (monosynaptic) or indirect (polysynaptic) (Varela et al., 2001). The connections between neural assemblies can be measured as the statistical dependence between two time-courses measured with functional MRI (fMRI) (Bolt et al., 2022). There is evidence that the transmission of information between distant neuronal assemblies is selectively mediated by different oscillatory frequencies that can be captured at the scalp-level with M/EEG techniques (Fries, 2015;Fries, 2005). Although there is some evidence for long-range synchronization being conveyed by high frequency oscillations (Varela et al., 2001;Arnulfo et al., 2020), typically, high-frequency oscillations modulate short-range local influences, while long-distance interactions are driven by slower rhythms (Canolty and Knight, 2010). This property falls within the scope of the standing wave model (Britannica, n.d.), which explains the generation of wave patterns in dynamical systems (Grassi et al., 2018). Recently, the standing wave-model has been suggested as the mechanism for zero-lag (i.e., instantaneous) functional activation (Bolt et al., 2022). According to this model, wave patterns arise from a direct association between the spatial structure and its temporal frequency of oscillation. A priori, all physical objects have a set of preferred frequencies (natural frequencies) at which the system tends to oscillate. Intuitively, one can refer to the oscillations of a violin string (Fig. 1), where the typical (resonant) frequencies of oscillations are determined notably by the length of the string; the longer the string, the lower the fundamental frequency. Similarly, ECoG studies have shown that, for distant electrodes in the brain, the correlation between the oscillatory activity is higher for low-frequency oscillations than for high-frequency oscillations (Raghavachari et al., 2006;Dickson et al., 2000;Canolty et al., 2007), and vice-versa. This suggests that, whereas low-frequency oscillations are better suited to modulate the amplitude of neural activity distributed over distant brain regions in longer temporal windows, high-frequency oscillations modulate activity over neighboring/local brain regions in shorter temporal windows (Canolty and Knight, 2010).
Novel tools derived from the field of graph signal processing (Grassi et al., 2018;Betzel et al., 2019) allow us to estimate the natural normal modes (spatial frequencies) of the connectome, coined connectome harmonics ( Fig. 1(c)). Connectome harmonics are expected to be the preferred frequencies at which the system tends to enter resonance (Atasoy et al., 2016), and when applied to EEG data, they have been proposed to capture the two main modes of neural information processing, namely localized and distributed, also referred as segregation and integration respectively (Wang et al., 2019;Glomb et al., 2020;Rué-Queralt et al., 2021). Low-frequency connectome harmonics correspond to patterns of activity associated with long-range connectivity, with a smooth change in value over distant nodes (Atasoy et al., 2016). In contrast, high-frequency connectome harmonics correspond to short-range connectivity patterns, i.e., only brain regions that are directly (and strongly) inter-connected show similar values. Consequently, it has been suggested that the dominance of low-frequency harmonics in EEG data reflects globally distributed (i.e., integration) mechanisms, whereas high-frequency harmonics capture localized activity (i.e., segregation mechanisms) (Wang et al., 2019;Glomb et al., 2020;Rué-Queralt et al., 2021). In functional MRI, the representation of BOLD signal activation in the connectome spectrum domain has been associated with the level of coupling, or dependency, of functional signals on the anatomical structure (Preti and Van De Ville, 2019). The well-spatially resolved, but slow integration time of the hemodynamic response function (Buxton et al., 2004) for an accurate mapping between the function and structure, but it is less suited to follow fast neural communication dynamics as M/EEG. When first discovered, connectome harmonics were shown to form the building blocks of well-known functional brain networks associated with both resting-state activity and different cognitive tasks (Atasoy et al., 2016). Recent work suggests that connectome harmonics can predict the specific temporal frequency range of a connectivity network in both fMRI (Atasoy et al., 2018) and MEG (Tewarie et al., 2019) recordings, indicating that the temporal synchronization between neural populations within specific connectome harmonics (i.e., graph frequencies) results in temporal-frequency-specific functional connectivity patterns. This link would imply that, at the whole brain-level, different temporal frequencies support information processing at different spatial The harmonic modes of a violin string refer to the different ways that a string can resonate when a force is applied. Each harmonic mode results in a different musical note. These modes are determined by the length of the string and its tension. The first harmonic mode, also known as the fundamental frequency, is the lowest frequency that a string can produce. This is the note that is played when the string is plucked or bowed in its natural state. The second harmonic mode is the next highest frequency that the string can produce, and so on. As the string length increases, the frequency of the harmonic modes decreases, resulting in a lower pitch. (b) Time. A circular graph is often used in stationary time series analysis because it provides a visual representation of the cyclical patterns in the data. In Fourier analysis, the harmonic modes of a circle refer to the different ways that a circular shape can oscillate or vibrate. (c) Connectome. In non-Euclidean domains, such as in graphs, the graph Fourier transforms has the same functionality, i.e., to represent the graph-signal into its main oscillatory modes. The Fourier or harmonic modes of the brain connectivity graph are known as "connectome harmonics".

scales.
Consequently, we aim at clarifying whether different temporal frequency bands of an EEG recording map into spatially distinct connectome harmonics during sensory perception. To test the validity of this hypothesis, we investigated well-known experimental paradigms, in which the neurophysiological processing of information in different frequency-bands has been studied invasively in animal models (Bastos et al., 2015). To this end, this work introduces the joint time-vertex connectome spectral analysis to EEG recordings during an auditory steady-state response (ASSR), which elicits evoked oscillatory activity time-locked to the stimulus onset (Grent-'t-Jong et al., 2021;Mancini et al., 2022a), and a visual grating task known to elicit an induced oscillatory response in several frequency bands (Hoogenboom et al., 2006;Mancini et al., 2022b) (Fig. 2(a)).
Based on the same underlying physical principles underlying the oscillating string ( Fig. 1), i.e., the wave equationwhich states that the frequency of resonance is inversely proportional to the length of the string-, we hypothesized that high-frequency oscillations should be associated with short-range functional connectivity (high-frequency harmonics, i.e., segregation mechanisms), whereas lower frequency oscillations should be associated with long-range connectivity patterns (low-frequency harmonics, i.e., integration mechanisms). It follows from this hypothesis, that the brain activity signal should have a compact representation in the joint time-vertex spectral domain. To validate this hypothesis, and the neurophysiological relevance of the joint time-vertex representation, we further investigated the correlation between the estimated patterns of localized and distributed neural activity in the temporal frequency domain and participants' behavioral performance during the visual task.

Recruitment and assessment of participants
The participants included in this study were the control subjects from a larger study from the 22q11DS Swiss Cohort (Schaer et al., 2009;Mancini et al., 2020;Bagautdinova et al., 2021). For this study subjects had to be aged between 17 and 30 years and should not have had any past or present neurological or psychiatric disease, use of psychotropic medications, psychopathology, learning difficulties, or premature birth. Further exclusion criteria specific to the task were the presence of auditory impairments documented by audiometry screening for the auditory task and a medical history of epilepsy or epileptic seizures for the visual task. A total of 32 participants met the inclusivity criteria in the present study, of which 30 in the visual study and 25 in the auditory study, with 23 subjects being included in both studies.
Written informed consent was obtained from participants and/or their parents. The study was approved by the cantonal ethics committee for research and conducted according to the Declaration of Helsinki.

Auditory steady state response paradigm
One hundred 40-Hz amplitude-modulated sounds (ripple tones: 1000 Hz carrier tones, duration 2 s) and 10 semi-randomly intermixed sounds without amplitude modulation (flat tones at 1000 Hz) were presented binaurally using intra-aural insert earphones (Etymotic Research, Elk Grove Village, Ill.). The interstimulus interval (ISI) was 2 s on average (1.5-2.5 s, uniform distribution). In order to ensure attentional engagement, participants were asked to detect the flat tones by button press and ignore the 40-Hz ripple tones that entrained the gamma-band responses (Glomb et al., 2020;Rué-Queralt et al., 2021).

Visual paradigm
The visual paradigm consisted of a centrally presented, circular sine wave grating. The circular grating drifted inwards toward the fixation point position and the speed of this contraction increased (velocity step at 2.2 deg/s) at a randomized time point between 750 and 3000 ms after stimulus onset (Preti and Van De Ville, 2019;Grent-'t-Jong et al., 2021;Mancini et al., 2022a). The experimental protocol consisted of 240 trials divided into three runs of 80 trials. Participants were instructed to press a button as soon as they noticed a speed increase. Stimulus offset was followed by a period of 1000 ms during which subjects were given visual feedback depending on their response. Before beginning the experiment, all the participants underwent a training session with one researcher to be sure that they understood the task. Behavioral measures were calculated as the percentage of correct answers out of the 240 trials and average reaction time. Participants are asked to report the change in speed of inward motion of the grating by button press. Right: 40 Hz Auditory steady state response (ASSR) paradigm composed by 100 ripple tones and 10 intermixed flat tones, each lasting 2 s. Participants were asked to recognize flat tones only responding by button press. (b) Electrical source imaging. Scalp patterns before, and after the sensory stimuli are presented (at 0 ms) and their source reconstruction at three time-points. (c) Connectome spectral analysis: gray matter areas defined by an anatomical parcellation define the nodes of the brain graph, and the estimated white matter tracts from diffusion MRI define the connectivity strength of the edges of the graph. From the constructed graph, we perform the graph Laplacian eigendecomposition and obtain a graph spectrum, consisting of a set of eigenvalues and a set of eigenvectors, the former ones termed graph spectrum, and the later termed graph or connectome harmonics. The source-reconstructed time-series signal represented in each brain area (top plot) or in each connectome harmonic (bottom plot).

EEG data acquisition and pre-processing
EEG data were continuously recorded with a sampling rate of 1000 Hz using a 256-electrodes HydroCel cap (Magstim-EGI) referenced to the vertex (Cz). The impedance was kept below 30 kΩ for all electrodes and below 10 kΩ for the reference and ground electrodes.
The preprocessing steps were performed by using the free academic software Cartool (Brunet et al., 2011)(https://sites.google.com/site/ca rtoolcommunity/home) as previously described (Mancini et al., 2022a(Mancini et al., , 2022bCantonas et al., 2021). In detail, the number of electrodes was reduced from 256 to 204 by eliminating the noisy electrode signals from the cheeks and neck. The data were band-pass filtered between 1 and 140 Hz using non-causal Butterworth filters and we applied additional notch filters at 50, 100, and 150 Hz. The periods of artifacts (e.g., muscle contraction) were marked manually and excluded from further analysis. Noisy channels were identified by means of visual inspection and excluded. Independent Component Analysis was used to remove eye movements (eye blinks and saccades) and ECG artifacts components using a Matlab script based on the EEGlab runica function EEGLAB (https://sccn.ucsd.edu/eeglab/) (Delorme and Makeig, 2004). The identified noisy channels were interpolated using a 3D spline interpolation (Perrin et al., 1989). Finally, we applied an instantaneous spatial filter implemented in Cartool that removes local outliers by spatially smoothing the maps without losing their topographical characteristics (Michel and Brunet, 2019) and data were recalculated to the common average reference.

MRI data acquisition and pre-processing
We acquired T 1 -weighted images with a 3-T Siemens Prisma at Campus Biotech in Geneva. The parameters for the acquisition of structural images for the T 1 -weighted MPRAGE sequence were as follows: TR = 2500 ms, TE = 3 ms, flip angle = 8 • , acquisition matrix = 256 × 256, 192 slices with field of view = 23.5 cm and voxel size = 0.9 × 0.9 × 1.1 mm, T 1 -weighted images underwent fully automated image processing with Freesurfer, version 6 (Fischl, 2012), comprising skull stripping, intensity normalization, reconstruction of the internal and external cortical surface, and parcellation of cortical and subcortical brain regions based on the Desikan-Killiany parcellation (Desikan et al., 2006). Then average measures of volume were extracted from those regions.

EEG analysis
For the visual task, only epochs with correct behavioral responses were considered for further EEG-analysis to ensure that the participants were paying attention to the task (Grent-'t-Jong et al., 2021;Mancini et al., 2022a). Time-frequency (TF) analysis was performed by using Morlet-wavelet transform (frequencies from 2 to 120 Hz, centered on steps of 2 Hz, with adapted resolution according to the FWHM (full width at half maximum) scheme) in MATLAB 2018b as described in previous papers (Mancini et al., 2022a(Mancini et al., , 2022b. Time epochs from − 1.5 to +1.5 s relative to the stimulus onset were averaged to event-related spectral perturbation (ERSPs) (Neuper and Pfurtscheller, 2001).
For source analysis, the inverse solution (IS) was computed using the academic software Cartool (Michel and Brunet, 2019) based on individual T1-weighted images preprocessed in Freesurfer. An approximate number of 5000 solution points were distributed in the individually segmented gray matter mask. We used the Locally Spherical Model with Anatomical Constraints (LSMAC) method for the lead field computation that was age-adjusted to reflect differences across age in skull conductivity and thickness (Michel and Brunet, 2019). A distributed linear inverse solution (Local AutoRegressive Average; LAURA) was employed to compute a transformation matrix from sensor level to inverse solution (Grave De Peralta Menendez et al., 2004).
The individual Desikan-Killiany parcellations obtained from Freesurfer were natively aligned on the brain of each individual and then used to label the 5000 solution points from the IS model in 84 ROIs covering cortical structures. Using the individual IS models, TF decomposition data from the surface were projected to the source space level, reduced as the power, normalized by the baseline period (− 1.5 to − 0.3 s ) and gathered as the median of the solution points in each ROI.

Consensus structural connectome
In order to formulate group predictions, we need to represent the brain activity signal of all participants on the same basis. For that purpose, and to obtain a reliable and robust connectivity estimate (Bolt et al., 2022), a unique group-consensus structural connectome was computed from a group of 70 healthy young adults (34 females), from an online dataset .
The mean age of the participants was 29.7 (range 18.5-59.2 years). All participants were scanned in a 3-Tesla MRI scanner (Trio, Siemens Medical, Germany) with a 32-channel head-coil undergoing diffusion and T1 weighted imaging.
The acquisition consisted in a diffusion spectrum imaging (DSI) sequence with 128 diffusion-weighted volumes and a single b0 volume (maximum b-value 8000 s/mm 2 ), using 2.2 2.2 3.0 mm voxel size. The DSI data reconstruction protocol is described in (Wedeen et al., 2008).
A magnetization-prepared rapid acquisition gradient echo (MPRAGE) sequence sensitive to white/gray-matter contrast (with 1 mm in-plane resolution and 1.2 mm slice thickness) was also applied, and Freesurfer (Delorme and Makeig, 2004) and Connectome Mapper 3 (Daducci et al., 2012;Tourbier et al., 2020) were used to segment the gray and white matter from the MPRAGE volume and to parcellate the gray-matter with the Desikan-Killiany parcellation (Desikan et al., 2006).
Deterministic streamline tractography on the reconstructed DSI data was performed to estimate the individual structural connectivity matrices, by seeding 32 streamlines per diffusion direction in each white matter voxel. The structural connectome was then defined as the number of fibers connecting each pair of brain regions given by the same parcellation used for the MPRAGE data. Within each brain region, the number of fibers found at the gray matter/white matter interface was summed.
From the connectomes of 70 healthy participants, a grouprepresentative consensus structural brain connectivity matrix was generated as in (Betzel et al., 2019). This method preserves the connection density of each single subject connectome, independently for intra-and interhemispheric connections. By doing so, this method allows more inter-hemispheric connections to be kept in the group estimate, in comparison to simple connectome average across subjects. The connection density of the consensus connectome is approximately 25%.

Connectome spectral analysis
The connectome spectral representation provides a different view of the data, just as PCA with all its components, but instead of this representation being data-driven, it is model-driven (by the connectome-spectrum). Connectome spectral analysis explains the brain activity signals in terms of their relationship with structural connectivity. The structural connectivity can be decomposed into Fourier modes, defining the so-called connectome harmonics (Atasoy et al., 2016). This is the graph-equivalent transform to the Fourier decomposition of temporal series (Atasoy et al., 2018).
Given the parcellated source reconstructed signal in time and frequency (s t,f ), we can obtain the contribution of each brain harmonic (ŝ t,f λ ) as: Where U r,λ indicates the weight of the λ-th connectome harmonic at the r-th brain region. The inverse transform, is very simple to obtain: And expresses brain activity at each time point as a weighted combination of coefficients that define the contribution of each harmonic.

Broadcasting index
Connectome harmonics are intrinsically ordered according to their spatial frequency in the connectivity graph (from low-to high-frequency harmonics). Originally inspired from the "structural-decoupling index" (Buzsáki and Draguhn, 2004), the broadcasting index (BI) is a temporally and time-frequency resolved real-valued index that indicates whether the overall-brain activity is polarized towards low-or high-frequency connectome harmonics, i.e., towards localized or distributed processing patterns.. The structural-decoupling-index describes how the functional activity time-series of a given brain region is "coupled" to the underlying connectivity structure of the brain. The BI is a temporally and time-frequency resolved real-valued index that indicates whether the overall-brain activity is polarized towards low-or high-frequency connectome harmonics, i.e., towards integration or segregation. Specifically, for each time-point we compute: L ) correspond to the proportion of the source reconstructed signal contained within high-(and low-) frequency connectome harmonics: where 1 [λ≥λT ] is an indicator function determining the connectome harmonics included in the norm-reduction. The dichotomization between high and low harmonics is based on the temporal statistics of the data, by selecting a threshold that splits the total energy of the signal into two equal halves (λ T ), for the entire data. See the original implementation (Rué-Queralt et al., 2021) for further details.
Positive values of BI indicate that the source reconstructed EEG signal is dominated by high-frequency connectome harmonics, thus by patterns composed of non-smooth terms more prone to present activity in localized brain regions. On the contrary, negative values of BI indicate the domination of low-frequency harmonics, thus by patterns composed of smooth terms encompassing a large number of brain regions (largescale networks).
The significance of the group's broadcasting index was assessed by a two-staged statistical test. In the first step, a non-parametric test against the null hypothesis of an average BI equal to zero. A null distribution was created by randomly permuting the individual's BI sign, at each time and temporal-frequency bin. A cluster-level correction was applied to the extracted p-values (Maris and Oostenveld, 2007).
A second test is then performed through Monte-Carlo simulations. Specifically, we designed a permutation test to generate a null model distribution of "no effect of structural connectivity on the broadcasting index". This null model was generated by computing the broadcasting index for n perm connectome surrogates. n perm was computed taking into account a significance level of α = 0.05, and the number of comparisons to be corrected for by the Bonferroni method, as follows: The number 50 was chosen as a trade-off to strike a balance between increasing the number of Monte Carlo samples for a more accurate pvalue estimation and preserving computational resources. Connectome surrogates consist of randomized graphs obtained using the Brain Connectivity Toolbox (Daducci et al., 2012) function null_model_und_sign. The harmonics from these graphs were obtained using the same approach as for the consensus connectome, and they are later referred to as surrogate harmonics in this manuscript.
The estimated broadcasting index was considered significant if it was within 5% of the extremes of the null model distribution. For this purpose, the uncorrected p-value was computed for each time-point as: Finally, we corrected the obtained p-values for multiple comparisons with the Benjamini-Hochberg FDR method (Benjamini and Hochberg, 1995), using the Python toolbox statsmodels.

Joint time-vertex spectral representation
The joint time-vertex spectral representation is a linear transform that consists of expressing brain activity (which is a function of location in the brain and time) in terms of (dual) frequency content (i.e. temporal frequencies and connectome harmonics). In this work, the time-spectral conversion is implemented via the Morlet-wavelet transform, to achieve temporally resolved time-spectra whereas the graph Fourier transform is used for the representation in the connectome graph. Given a matrix containing in each element S r,t the activity signal for the brain region r at time t, the joint time-vertex representation is obtained as: Ŝ is joint time-vertex Fourier transform of S, where U is linear operator corresponding to the graph Fourier transform (see Connectome Spectral Analysis) and W is the linear operator that maps each brain activity at each time-point to a vector of Morlet-wavelet coefficients.

Compactness of signal
The compactness of a signal is an estimator of its sparsity i.e., how many zero-coefficients it has. In this work, we use the L1/L2 ratio to measure compactness of a signal: r = ‖s‖ 1 /‖s‖ 2 where the L2 norm (square root of the sum of the squared value of the elements) in the denominator is used to compare the L1 norm (sum of absolute value of the elements) among signals of different energy. If we compare two signals with the same total magnitude (here estimated with the L2 norm), the sparser signal (signal with more elements zeroed, or close to zero) will have a smaller L1 norm. A low value of the L1/L2 ratio indicates that the signal has many low-valued coefficients (compactness) whereas a high value indicates the presence of many high-valued coefficients.
When a transformation represents the degrees of freedom in the data using a smaller number of non-zero coefficients (in a more concise manner), it indicates a strong parametrization. This idea aligns with the fundamental concept behind PCA analysis, which aims to capture a significant portion of data variance with the fewest components possible, but without relying on statistical properties of the data.

Joint time-vertex spectral representation of brain activity
We transformed the source-reconstructed EEG signal, at each time point, jointly by the Morlet-wavelet transform (time domain) and the connectome Fourier transform (brain connectivity domain), yielding a weighted combination of active temporal frequencies and connectome harmonics (see Methods: Joint time-vertex spectral representation).
Our first observation was that the joint time-vertex spectral representation is the most compact, compared to just the vertex-or the timespectral representation: the compactness of each signal was assessed by the ratio (median 1075.30, vs. 1525.25 for Morlet Wavelet representation, Wilcoxon signed rank test p-value<0.001). This result indicates that the pattern of brain activation can be understood in terms of the involvement of very few and specific brain networks (compact, specific representation), which is more informative than the distributed activation pattern observed in the brain region representation. We isolated the brain networks underpinning the dynamic responses to the different stimuli observed in Fig. 3(b,b'). These brain networks captured (individually) more than 30% of the variance of a given oscillatory rhythm (Fig. 4).
In the visual task, at around 0.2 s after stimulus presentation, there is a peak in high gamma activity (58-68 Hz) corresponding to activation of the primary visual cortex and surrounding areas. At about the same time there is a trough in the alpha (9-12 Hz) and beta (15-30 Hz) bands reflecting desynchronization (Fig. 3(c)), corresponding to a wider fronto-parietal network (Fig. 3(d)).
In the auditory task, oscillatory activation is time-locked to the stimulus onset and corresponds to the activation of temporal regions and the anterior cingulate cortex (ACC) for gamma oscillations (38-42 Hz) and a wider fronto-parieto-temporal network for theta oscillations (4-8 Hz) (Fig. 5(d)).

Frequency bands as channels for broadcasting brain activity
Next, we computed the broadcasting index (BI, see Methods) of brain activity for each time point and temporal frequency, following the methodology we introduced previously in (Rué-Queralt et al., 2021). Briefly, the BI summarizes, in a single value, the connectome harmonic composition of the source-reconstructed EEG. A positive BI value indicates that the EEG signal pattern has a bias toward connectome harmonics with large eigenvalue (high graph-frequency), which defines an important contribution of the brain networks supporting local neural processes. On the contrary, a negative BI value indicates a low graph frequency dominance reflecting distributed coordinated neural activity. Figure 5(a,d) show how, in both visual and auditory tasks, there is a dominance of smooth connectome harmonics in the alpha and theta frequency bands (right after stimulus presentation). On the other hand, the more localized neural activity appears to be specific to higher temporal frequencies, namely high gamma (60-70 Hz) in the visual task and low gamma (40 Hz) in the auditory task.  Figure 5(b,e) illustrates how BI varies as a function of temporal frequency band. Specifically, it shows the conditional probability of a BI value given the activation pattern at a frequency in Hz. Brain activation patterns at low temporal frequencies (e.g., up to 30 Hz) occur more likely in the negative BI value range, indicating the dominance of distributed neural processes mediated through long-distance connectivity. At the higher end of the temporal frequency spectrum, activity tends to have a much smaller negative BI magnitude, indicating that this frequency range codes primarily local neural processes mediated through shorter range connectivity.
Based on our results showing the dominance of specific temporal frequency bands in the neural response toward different stimuli (Fig. 2), we summarized the respective BI values for theta, alpha, beta and gamma (low and high) bands. Figure 5(g) shows the distribution of BI values across these bands, indicating similar BI specialization in response to visual and auditory tasks. The pattern of significant differences among the BI values in different temporal frequency bands suggest the existence of two main frequency-band clusters: the first consists of the theta and alpha bands, which are associated with distributed neural activity (large and negative BI values), whereas the second group consists of beta, low-gamma and high-gamma (small and negative BI values), and are linked to more localized neural activity (Wilcoxon ranksum's test p-value<0.05, corrected for multiple comparisons with FDR, see Fig. 5(h,i)).

Broadcasting index predicts behavioral performance
During the visual stimulus task, we recorded the percentage of accurate responses for each participant and reported their distribution in Fig. 6(a). We tested whether behavioral responses were related to the BI. For each time point and temporal frequency, we calculated the Spearman's correlation coefficient between the BI value and the mean behavioral score (detection rates) of all participants (see Fig. 6(b)). The results show an initial distributed neural activity in the alpha band during the first 300 ms, followed by an increase in localized processing (towards less negative BI) in the lower gamma range (Spearman's correlation test compared against null distribution of shuffled participants' score, p-value<0.05, corrected for multiple comparisons with Bonferroni method). This contrasts with the more generalized pattern of correlation of the performances with signal power, which partially overlaps (see Fig. 6(d)). Suggesting that higher performance is associated with a predominance of distributed activity in the lower temporal frequencies and of localized processing in the higher ones.

Discussion
Recent efforts to use anatomical connectivity priors to estimate neural functional connectivity have shown improved estimation for electrophysiological signals  and fMRI data (Crimi et al., 2021). In this work, we used anatomical connectivity to estimate the frequency-resolved dynamics of different modes of information processing in the brain and validate it with existing literature using intra cortical recordings. Specifically, our method estimates the different modes of information processing by describing the source-reconstructed EEG signal not only in terms of connectome harmonics (Atasoy et al., 2016), i.e., the natural normal modes (spatial frequencies) of the brain connectome ( Fig. 2(c)), but also with respect of its temporal frequency content. This is achieved by using a joint time-vertex spectral representation. In this representation, the pattern of brain activity is understood as a weighted combination of the active connectome harmonics at different temporal frequencies.
We have applied this joint time-vertex spectral representation of brain activity on two different datasets, an auditory steady state response (ASSR), which elicits evoked oscillatory activity time-locked to the stimulus onset (Mancini et al., 2022a), and a visual grating task known to elicit induced oscillatory response in several frequency bands (see Fig. 2). We found that the joint spectral representation has a highly structured pattern (compact and ordered), showing activation for specific connectome harmonics and temporal frequency pairs (i.e., for specific spatial and temporal scales). This result indicates that the joint time-vertex spectral representation captures the major axes of variability in the data. These axes of variability, or degrees of freedom, are graph waves and temporal oscillations and they are tightly linked.
Importantly, similarly structured time-vertex spectra are also encountered in physical phenomena governed by the wave equation (Grassi et al., 2018). Our results suggest that brain electrical activity follows a similar principle, at least in the situation where there is significant functional coupling mediated through structural connectivity, as in the setting that has been studied here. Thus, we add evidence to the notion that brain spatio-temporal oscillatory patterns are profoundly shaped by the underlying topology of the connectome. Although this result might be expected in part from the literature showing that functional connectivity networks can be sparsely reconstructed by connectome harmonics (Wang et al., 2019;Glomb et al., 2020;Rué-Queralt et al., 2021) or through computational models (Lakatos et al., 2005), to our knowledge, this is the first evidence for compactness, in terms of ratio, in the joint time-vertex spectral domain.
Using the joint-time vertex spectral representation, we isolated the connectome harmonics that underpin the dynamical response to the different tasks (Fig. 5). Looking at the main networks involved in the perceptual processes, we found consistent patterns in the two experiments. The peak gamma activity corresponded to the activation of relatively high frequency connectome harmonics that captured the primary sensory cortex and the surrounding areas, i.e., the areas of the early visual system for the visual task, and the temporal regions and the anterior cingulate cortex for the auditory task. In contrast, at lower temporal frequencies, we found dominance in smoother connectome harmonics of broader spatial distribution typically known to be involved in more complex cognitive processes, namely a fronto-parietal network in the alpha and beta-bands for the visual task, and a fronto-parietotemporal network in theta for the auditory task. These results also agree with the current literature (Pagnotta et al., 2022;Costa-Faidella et al., 2017) and validate the representation of brain activity in the connectome spectral domain.
We used the joint time-vertex spectral representation, to estimate the broadcasting index (BI) of the source-reconstructed EEG data. The BI describes, at a given point in time and temporal frequency range, the connectome harmonic composition of the source-reconstructed EEG signal, indicating the dominance of localized or distributed neural processing. Inspired by recent developments in the field of graph signal processing (Grassi et al., 2018), our work generalizes the concept of BI dynamics (Rué-Queralt et al., 2021) to the temporal frequency domain. A negative BI indicates that low-frequency harmonics are more active than high-frequency harmonics, suggesting a predominance of widespread coordinated neural activity, which may refer to the well-known integration mechanisms of the network neuroscience literature (Fries, 2015). On the other hand, positive values suggest that local coordinated activity is dominant, which may be related to functional segregation. Our results show a very strong relationship between temporal frequencies and BI, further supporting the idea that the frequency of synchronization between brain regions -which can be understood according to the communication through coherence concept (Fries, 2015;Fries, 2005) as a form of neural communication -is intimately linked to the underlying connectivity structure (Fig. 5).
For both tasks, the BI had its largest values at high and low gamma, indicating that these temporal frequency bands tend toward local synchronization, i.e., functional segregation. Theta and alpha show more negative values, indicating a tendency toward large scale activity, supporting functional integration. While beta band activity has, in general, a balanced composition of connectome harmonics for the visual task, it is associated with slightly more low frequency harmonics in the ASSR task. This close relationship between information processing modes and oscillation frequency is again tightly linked with the wave equation model, as previous studies had suggested (Roberts et al., 2019;Huntenburg et al., 2018;Breakspear et al., 2006;Raut et al., 2020). Using a computational model of brain activity propagated through the structural connectome, the authors of (Roberts et al., 2019) suggest that brain activity could be composed of wave-like patterns that are reconfigured over time. Recent studies further suggest a functional specialization of these brain waves (Huntenburg et al., 2018;Raut et al., 2020), as well as a dysfunctional involvement in the development of epileptic seizures (Breakspear et al., 2006).
As said, there are studies suggesting that the transmission of information between distant neuronal populations is mediated by different oscillatory frequencies (Fries, 2015;Fries, 2005). At the network level, this would translate into networks with long average connectivity (low-frequency harmonics in charge of integration mechanisms) oscillating at low temporal frequency (theta and alpha), and networks with short average connectivity (high-frequency harmonics in charge of segregation mechanisms) preferentially oscillating at high temporal frequency (gamma). Computational neuroscience models support this hypothesis (Raj et al., 2020), which is in agreement with our results. Overall, our findings suggest that the application of connectome spectral analysis could yield new insights into this theory by identifying the (temporal) spectral composition of each (spatial) harmonic under different experimental conditions or simulations. It should be noted that high-frequency oscillations have been also implicated in long-range synchronization (Arnulfo et al., 2020). However, synchronization in the gamma range is strongly dependent on anatomical cortico-cortical connections (Vezoli et al., 2021) and preferentially conveys the signal in a feedforward direction. Therefore, gamma rhythm originates locally and propagates following predefined anatomical routes in a bottom-up fashion, making it less suitable to support large-scale cortico-subcortical polysynaptic networks implicated in integration mechanisms.
Finally, we investigated whether the broadcasting index (BI) could predict behavior. We found that the performance in the visual task is associated with greater localized activity in the gamma range and greater large-scale activity in the alpha band range (Fig. 6). These results suggest that the brain state resides in a well-balanced sequence of (c) Probability densities of the variance of task accuracy explained by the broadcasting index (blue) and signal power across brain regions. The distribution is defined across all bins showing significant correlation in time and frequency. (d) Time-frequency plot of the bins in which Spearman's correlation between task accuracy and BI (blue) and power (pink) are significant. Those time-frequency bins that are significant for both BI and power are colored in black. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) specialized activity at the level of the visual cortex and integration in the prefrontal and associative regions. Alpha desynchronization is thought to mediate top-down information processing, reflecting attentional processes toward visual stimuli (Rihs et al., 2009;Romei et al., 2008), and thus our findings may reflect the top-down propagation of alpha oscillations across different brain regions to achieve maximal integration. On the other hand, it has been proposed that gamma-band oscillations transmit the propagation of sensory stimuli in the respective sensory cortices (Fries, 2015). Therefore, greater segregation in the gamma band range may indicate greater efficiency in local processing at the level of the visual cortices.
In summary, this study presents a novel approach to explore new hypotheses on brain communication mechanisms linked to the brain connectivity structure. We analyzed the complex interplay between distributed versus localized neural activity in the temporal frequency domain, in the specific cases of visual and auditory perception. Our findings suggest that different communication mechanisms take place in different temporal frequency channels broadcasting information at different spatial scales, and that they predict the performance in a behavioral task. Specifically, in line with current knowledge from invasive studies (Bastos et al., 2015;Vezoli et al., 2021), we found that high frequency oscillations are implicated in localized, i.e., segregation, neural mechanisms while low frequency oscillations in the integration of sensory information. Thus, we have demonstrated a principled approach to study spatio-temporal coupling of neural activity which may contribute to better understand different modes of information processing not only in health but also neuropsychiatric disorders (Lord et al., 2017) or other brain disorders such as epilepsy (Bastos et al., 2015).

Limitations and future research
The study of brain activity using connectome harmonics has some limitations related to the fact that diffusion MRI-based tractography contains artifactual connections. However, this is not a crucial limitation since the most important variability in connectome harmonics comes from the local connectivity of the brain (Naze et al., 2021).
Another limitation is the use of spatial smoothing priors in the source reconstruction algorithm, as this imposes a small but consistent shift of energy towards the low frequency connectomic harmonics (Rué-Queralt et al., 2021). Our method of estimating the broadcasting index aims at being agnostic to this bias by choosing a threshold on the connectome spectrum that divides the energy of the integration and segregation into two equal parts. In this way, if the energy is shifted to one end of the spectrum, the threshold is moved accordingly.

Conclusions
In summary, this study uses a novel method to explore new hypotheses on brain communication mechanisms linked to the brain connectivity structure. We analyzed the complex interplay between distributed and localized neural processes as a function of temporal frequency in the specific cases of visual and auditory perception. Our findings suggest that different communication mechanisms take place in different frequency channels broadcasting information at different spatial scales, and that they predict the performance in a behavioral task. Specifically, in line with current knowledge from invasive studies (Bastos et al., 2015;Vezoli et al., 2021), we found that high temporal frequency oscillations are implicated in segregation while low temporal frequency oscillations in the integration of sensory information.

Data and code availability
Data and code will be made available upon reasonable request by contacting the corresponding authors.

Declaration of Competing Interest
Authors declare that they have no conflict of interest..

Data availability
Data will be made available on request.