Cross-Frequency Multilayer Network Analysis with Bispectrum-based Functional Connectivity: A Study of Alzheimer’s Disease

Alzheimer’s disease (AD) is a neurodegenerative disorder known to affect functional connectivity (FC) across many brain regions. Linear FC measures have been applied to study the differences in AD by splitting neurophysiological signals, such as electroencephalography (EEG) recordings, into discrete frequency bands and analysing them in isolation from each other. We address this limitation by quantifying cross-frequency FC in addition to the traditional within-band approach. Cross-bispectrum, a higher-order spectral analysis approach, is used to measure the nonlinear FC and is compared with the cross-spectrum, which only measures the linear FC within bands. This work reports the first use of cross-bispectrum to reconstruct a cross-frequency FC network where each frequency band is treated as a layer in a multilayer network with both inter- and intra-layer edges. An increase of within-band FC in AD cases is observed in low-frequency bands using both methods. Cross-bispectrum also detects multiple cross-frequency differences, mainly increased FC in AD cases in δ-θ coupling. Overall, increased strength of low-frequency coupling and decreased level of high-frequency coupling is observed in AD cases in comparison to healthy controls (HC). Vulnerability analysis reveals that the integration and segregation properties of networks are enabled by different frequency couplings in AD networks compared to HCs. Moreover, the segregation property depends less on the γ interactions in AD cases, suggesting a shift from high to low-frequency connectivity. The results highlight the importance of studying nonlinearity and including cross-frequency FC in characterising AD. Moreover, the results demonstrate the advantages and limitations of using cross-bispectrum to reconstruct FC networks.


I. INTRODUCTION
Alzheimer's disease (AD), the most common form of dementia, causes early degradation of neural circuits leading to cell death and synaptic loss [1]. Studies have shown that AD affects distributed brain networks, alters functional connectivity, which can lead to disconnection syndrome and disrupts information processing across multiple scales [2][3][4][5].
Electroencephalography (EEG) is a common method to study AD. Currently EEG is not commonly used for the diagnosis of AD, but many studies demonstrate a high potential in developing EEG-based biomarkers and diagnostic tools. The main EEG characteristics associated with AD are the slowing of signals and altered synchronisation [3,[5][6][7][8]. Slowing of EEG in AD was observed as increased power in δ and θ frequency bands and decreased power in α and β frequency bands [5,6]. Similarly, AD shows increased synchronisation within lowfrequency bands (<12 Hz), decreased synchronisation in high-frequency bands, and is associated with the altered functional connectivity, especially the long-distance cortical connections [9]. However, these characteristics are typically measured at a single channel or between channel pairs. In contrast, network-based methods analyse multiple channels and reveal additional characteristics of * Correspondence to: fei.he@coventry.ac.uk AD, namely, reduced integration of information [10,11] and loss of small-worldness [12]. However, these characteristics are often analysed only within specific frequency bands.
This study aims to extend the functional connectivity beyond within-frequency coupling (WFC) taking the cross-frequency coupling (CFC) [13] into account. WFC networks of AD were analysed previously by using coherence (linear) [14] and wavelet coherence (nonlinear) [15]. Only one CFC measure, i.e. phase synchronisation index (PSI), had been used for the graph analysis of CFC networks in AD [16]. This work extended the findings of reduced integration and loss of small-worldness to CFC multilayer networks. However, it does not consider the roles of different frequency components in the networks. The multilayer-network framework had been used previously for brain network analysis. Loss of interfrequency hubs in AD was reported using MEG multiplex networks [17,18], where the inter-layer edges are inserted with fixed weight only between the same nodes across layers. Alterations in multilayer network hubs have been reported in multimodal networks in AD [19] and fMRI frequency-band networks in schizophrenia [20]. Multilayer networks integrating WFC and CFC have been used to analyse MEG data from healthy [21] and schizophrenic subjects [22]. Tewarie et al. [21] show that frequencyband network layers interact via CFC, share a certain amount of common structure and operate at the edge of independence and inter-dependence. However, these For each pair of EEG electrodes, a cross-bispectrum is estimated. The frequency bands coupling edge weights are given by the maximum value within the respective CBS window, e.g. δ-δ (red). Note that CBS estimates are directed, e.g. δ-γ = γ-δ (both in grey). Thus from each CBS, 25 edges are inferred. (C) Using the edge weights inferred from CBS, a multilayer network is constructed with layers representing the frequency bands of EEG. Such a network has both intra-layer edges and inter-layer edges, representing the within-frequency coupling and cross-frequency coupling, respectively. studies analyse the layer-relationships mainly as correlation of their adjacency matrices or as differences in global average coupling strength.
Bispectrum is a higher-order spectral analysis and quantifies quadratic coupling between two frequency components and their algebraic sum [23].
It has been shown to detect amplitude-amplitude and phaseamplitude CFC in addition to the phase-phase coupling [24,25]. The bispectral coupling also indicates an increase in non-Gaussianity [26]. Features derived from bispectrum were proposed as biomarkers of epilepsy [27], Parkinson's disease [28] and AD [26,29,30]. Most of these studies compute (cross-) bispectra of only a few channels or pairs of channels. Although, a few studies used bispectrum to compute global networks, these analyses do not use graph-theory. Instead, each node is analysed in isolation [29] or single-channel bispectra are averaged across nodes to derive certain global properties [30]. In contrast, this study computes crossbispectra between all pairs of EEG channels to estimate the widely distributed FC brain networks and perform graph-theoretical analysis.
In this work, the cross-bispectrum (CBS) estimates of functional connectivity are computed. We would like to investigate the contribution of nonlinear WFC and CFC in differentiating between Alzheimer's disease (AD) and healthy controls (HC) in comparison to the equivalent linear WFC measured with cross-spectrum (CS) (Fig. 1). To the best of our knowledge, this is the first attempt to use CBS to reconstruct FC networks in patients with AD. Moreover, we report a novel multilayer-network-based approach to elucidate the roles of the traditional EEG frequency bands and their CFC in the sensor-level EEG networks of HC and AD.

II. DATA
EEG recordings were collected from 20 AD patients and 20 healthy control participants (HC) younger than 70 years. A detailed description of the experimental design and confirmation of the diagnosis is provided in [31]. All AD participants were recruited in the Sheffield Teaching Hospital memory clinic. AD participants were diagnosed between 1 month and 2 years before data collection, and all of them were in the mild to moderate stage of the disease at the time of recording with the average Mini Mental State Examination score of 20.1 (sd = 4). Age and gender-matched HC participants with normal neuropsychological tests and structural MRI scans were recruited.
EEG was acquired using an XLTEK 128-channel headbox, Ag/AgCL electrodes with a sampling frequency of 2 kHz using a modified 10-10 overlapping a 10-20 international electrode placement system with a referential montage with a linked earlobe reference. The recordings lasted 30 minutes, during which the participants were instructed to rest and not to think about anything specific. Within each recording, there were two-minutelong epochs during which the participants had their eyes closed (alternating with equal duration eyes-open epochs, not used in this work).
All the recordings were reviewed by an experienced neurophysiologist on the XLTEK review station with time-locked video recordings (Optima Medical LTD). For each participant, three 12-second-long artefact-free epochs were isolated. Finally, the following 23 bipolar channels were created:  A. EEG pre-processing EEG signals were confirmed to be artifact-free thus no additional artifact removal was undertaken. The signals were band-pass filtered to be between 0.1 and 100 Hz using a zero-phase 5 th order Butterworth filter; 50 Hz relating to the power line noise was removed using a zero-phase 4 th order Butterworth stop-band filter, and the data were down sampled to 250 Hz using an order 8 Chebyshev type I filter. Finally, the signals were normalised (to zero mean and unit standard deviation).

A. Cross-spectrum and cross-bispectrum
The spectrum of signal S X is calculated via smoothed periodogram. Fast Fourier Transform (FFT) is used to estimate the periodogram with Daniell smoothers. The periodogram is computed over 256 frequency bins (0.98 Hz bandwidth). CS at frequency f is then computed as: CS XY (f ) = S X (f ) · S Y (f ). An absolute value of CS is calculated. A direct FFT-based method is used to estimate absolute CBS: where · denotes averaging, X T (f ) is a Fourier Transform of signal X over an interval T and Y * is the complex conjugate. 256-point FFT is used. CBS is computed over 1-second-long segments with 50% overlap over whole frequency range (0.5 -100 Hz). The estimated CBS is then smoothed in the frequency domain using Rao-Gabr window (size 5). CS and CBS were computed for all pairs of EEG channels. Five frequency bands b are considered: The connectivity between channels X and Y and frequency bands b X and b Y is computed as: for CS and CBS respectively. This resulted in five WFC (CS and CBS) and 20 CFC (CBS only) measures for per channel pair. It is of note that the CBS is directed. In order to ensure statistical significance of the estimated connectivity, surrogate thresholding was used [32]. For each pair of channels, 200 surrogate signals were generated using the Fourier transform surrogate which scrambles the phase of the signal, and their CS and CBS are computed. The 95% confidence interval of surrogate values is computed and used as a threshold. Coupling values below the threshold are set to zero. We chose this approach to ensure the reliability of estimated brain networks. In contrast, Wang et al. [29] used no such thresholding when analysing bicoherence coupling. Alternative approaches exist in the literature. Chella et al. [33] take advantage of the asymmetric property of CBS to ensure robustness against mixing artifacts.
We obtain a set of connectivity matrices for each EEG recording, i.e. N × N matrices (N = 23). For CS and CBS, there are five and 25 connectivity matrices, respectively. A global (averaged per subject) connectivity is computed for each 23 × 23 matrix and compared between groups using a Mann-Whitney test.

B. Network measures
To identify the important channels in the network, we compute a coupling specific node strength (NS) [34] for each channel i and different types of frequency couplings c, where Π(i, c) are the nodes connected to channel i via edge type c and w ij is the edge weight, i.e. CS or CBS connectivity given by ijth entry of the N ×N connectivity matrix. This measure is computed for both CS and CBS, resulting in 5 (5 frequency bands) and 25 (5×5 frequency bands) values per channel, respectively. In order to analyse the importance of the different frequency couplings in terms of the whole brain network, we represent them as a multilayer network. In this network, nodes are located within layers representing the different frequency bands, WFC represents the edges between nodes within a single layer, i.e. intra-layer, and CFC represents the edges between nodes located in different layers, i.e. inter-layer. In this paper, the CS networks are not analysed as multilayer networks, since such networks would have no inter-layer edges and thus would not be comparable directly with the CBS networks. The following measures are computed only for CBS networks. We obtain networks with 23 nodes that are replicated over 5 layers (L ∈ [δ, θ, α, β, γ]), resulting effectively in 115 nodes. There are 5 types of intra-layer edges such as δ-δ, and 20 types of inter-layer edges, such as δ-θ or θ-δ.
We measure the importance of each type of frequency coupling within the multilayer network by measuring the contribution of each edge to enable efficient passing of information through the network. For this purpose, we define coupling betweenness centrality (CBW) based on an adjusted version of edge betweenness (BW) [35]: where E is the total number of edges of coupling type c and BW (e) is edge betweenness centrality given by: where g ij is the number of shortest paths (least sum of 1/w ij ) between nodes i and j, and g ij (e) is the number of those paths that go through edge e. CBW quantifies the contribution of each coupling type to the information integration [36], i.e. the amount of information flow through edges. CBW assumes that the essential processes within the network occur along the shortest paths. However, there might be alternative paths with only minor length differences, which CBW ignores. In case of a disruption of the network structure, these alternative paths might enable recovery of function with negligible differences. We quantify this as the vulnerability of the network to the removal of one type of frequency coupling. The vulnerability is measured in two ways, namely the loss of ability to integrate information [37] and the loss of segregation.
The integration property of network G, i.e. the ability of a network to communicate information globally, is approximated with global efficiency (gE) given by: where N is the number of nodes in network G and d(ij) is the shortest path length between nodes i and j. gE is related to CBW. CBW measures the information flow on the more detailed edge-level while gE takes the nodelevel perspective. The segregation property of network G, i.e. the presence of densely connected clusters and sparse connections between them, is approximated with a local efficiency given by: where G i is the neighbourhood of node i, i.e. subgraph of nodes directly connected to i, without node i itself.
In order to measure the vulnerability of the network and its dependence on different types of frequency coupling, gE and lE are computed for the full network. The two measures are then re-computed on a perturbed networks where one type of frequency coupling (i.e. set of edges) was removed. The change in gE and lE give the global and local vulnerability measures gV (G c ) = 1 − (gE(G c )/gE(G)) and lV = 1 − (lE(G c )/lE(G)), where G is the full network and G c is the perturbed network.

C. Network thresholding and statistical analysis
In order to filter out the unimportant edges that might be the result of a spurious coupling, the weighted multilayer networks are thresholded through relative quantilebased thresholding. Given a quantile Q, all edges with weight lower than Q are removed from the network. There are large differences between the weights of each frequency coupling type (e.g. mean of γ-β = 1.627 compared to mean of α-α = 8.975); thus, a separate threshold Q is used. As a result, the networks retain Q% of the strongest edges. To ensure that the observed differences between the networks are not due to the choice of threshold Q, all of the network measures are computed over 10 of threshold values (Q ∈ [0.5, 0.95] in increments of 0.05) and only significant differences observed over at least five thresholds are declared significant. The reported plots and numerical results are obtained with Q = 0.7 as this is closest to the median of the threshold values. Alternatively, Q = 0.75 could be used. However, the effect of choice of this threshold should be minimal as the significant differences were observed over multiple thresholds. We report the number of consecutive thresholds where significant differences were observed. All p-values are corrected using the Benjamini-Hochberg false discovery rate method [38].
As NS compares the different frequency coupling types in isolation, the weighted networks are used in this analysis. Furthermore, we convert the CBS network from directed to undirected by taking the mean weight for each pair of directed edges, thus collapsing them into a single edge. NS is log-transformed to reduce skewness. As NS is not normally distributed, a non-parametric Mann-Whitney U test is used. NS is computed for both CS and CBS networks and frequency bands that show a significant difference in the global averaged coupling strength.
The multilayer graph measures aim to analyse the roles of the frequency coupling types in the properties of the global network. Since differences were observed between the weights of frequency coupling types, we convert the weighted networks into unweighted networks by setting the weight of all edges to 1, thus removing any effect of the observed differences. CBW , gV and lV are not normally distributed, thus, the Mann-Whitney U test is used to compare between groups. Furthermore, CBW and gV were log-transformed to reduce skewness.
In order to minimise the effect of the selected artifactfree epochs, we average all of the compared measures described above over all three epochs for each subject. This is to ensure that the observed differences are independent of the epoch selection and can be observed robustly in an arbitrary time window during the EEG recording.

IV. RESULTS AND DISCUSSION
A. Connectivity matrices and average connectivity Averaged connectivity matrices ( Fig. 2 and 3) indicate that both methods seem to detect differences in the topology of functional connectivity networks. Using CS, significant differences in average connectivity were found in δ and θ where AD cases have increased connectivity. Using CBS, differences can be observed in multiple frequency bands and their couplings. Increased global connectivity is observed in AD cases in θ WFC and θ-δ CFC. In contrast, decreased global connectivity in AD cases is found in α WFC and α-β, α-γ and β-γ CFC. Overall, AD cases show increased connectivity of low-frequency components and their CFC interactions and decreased connectivity in high-frequency components These findings are consistent with the literature in reporting an increased activity in δ and θ in AD [3,5]. An increase in δ WFC and low-frequency CFC in AD was also reported using bicoherence [29]. Similarly, Maturana et al. [30] report increased bispectral power in AD in δ and θ and a decrease in α, β 1 and β 2 . They also report lower bispectral entropy in δ and θ suggesting less frequency components interact with these frequency bands. Furthermore, a decrease in γ-θ phase-amplitude coupling was seen to be a sign of progression from mild cognitive impairment to AD [39]. In contrast, Cai et al. [16] report the opposite differences in the same WFC and CFC using PSI, i.e. decrease in δ and θ.
Moreover, the visible distortion of structure within multiple frequency bands detected by both CS and CBS suggests connections to the disconnection syndrome and disturbed information processing in AD.

B. Coupling-wise node strength
In order to statistically test the differences in connectivity measured by both CS and CBS and to localise the brain regions which show the most pronounced differences between AD and HC, NS is measured for each channel and coupling type separately. Due to the limits, we do not report the numerical results of these tests. Alternatively, we show the results in Figs. 2 and 3 for CS and CBS, respectively.
The differences in WFC detected by CS and CBS ( Fig.  4 and diagonal elements in Fig. 5) are generally similar. Both methods show increased θ NS in AD cases across most channels. Both CS and CBS show decreased α in parietal and occipital channels. In addition, CBS shows differences in several channels that CS failed to detect, such as decreased γ in AD cases in frontal and temporal channels. On the other hand, CBS fails to capture the increased NS in δ in AD cases that can be seen in CS. These differences showcase the importance of assessing both linear and nonlinear coupling in understanding the variations in AD brain networks.
FIG. 6: Importance of each type of frequency coupling of HC (blue) and AD (orange) measured by edge betweenness. Significant differences (p ≤ 0.05) observed in at least five thresholded networks are encoded by asterisks.
Multiple differences in CFC (off-diagonal elements in Fig. 5) were detected, highlighting the need to analyse the relationships between frequency components both in healthy and AD brain networks. AD cases show a widespread increase in δ-θ and δ-α in frontal channels and their inverse. This is in contrast with Wang et al.
FIG. 7: Global vulnerability of HC (blue) and AD (orange). Significant differences (p ≤ 0.05) observed in at least five thresholded networks are encoded by asterisks. [29] where an increase of δ-θ only in frontal channels is reported. They also report an increase in midline parietaloccipital θ-γ that we did not detect. Furthermore, we observed a parietal-occipital and temporal decrease in αβ and frontal decrease in α-γ and β-γ in AD cases. Cai et al. [16] report comparable differences using PSI, TABLE I: Differences in coupling betweenness centrality (CBW), global vulnerability (gV) and local vulnerability (lV).
Comparisons were done using Mann-Whitney U test. The results are reported in the following order: test statistic, p-value and number of thresholds where significant differences were observed. The input frequency is on the vertical facets, and the output frequency is on the horizontal. Significant differences that were observed over at least 5 thresholds are in bold. but in contrast to our results, they show mainly decreased NS in AD cases. This might be because CBS is influenced by the amplitude, while PSI is a pure phase coupling measure. Fraga et al. [40] report an increase of the δ-θ and δ-β amplitude-amplitude CFC in AD cases which is similar to our results. This suggests that CBS indeed measures some mixture of CFC types [24], since our results are partially in line with phase-phase CFC and amplitude-amplitude studies.

C. Multilayer network analysis
In order to elucidate the roles of the frequency bands and their coupling, both WFC and CFC, we analyse the CBS networks as multilayer networks with 5 layers representing the traditional frequency bands of EEG.
First, CBW is used to assess the importance of each type of coupling for both local and global communication in the network. Numerical results of comparing CBW are reported in Table I and visualised in Fig. 6. AD cases show a decreased importance of α-γ, β-γ, γ-θ and γ-α. This finding is likely linked to the characteristic slowing down of signals in AD, i.e. decrease of highfrequency power [5,6]. On the other hand, we observe an increased importance of α-θ, β-θ, β-α and γ-γ in AD cases. Interestingly, the importance of γ-γ is also the lowest overall.
Then, gV is used to assess the vulnerability of information integration of the network to the removal of a coupling type. Numerical results of comparing gV are reported in Table I and visualised in Fig. 7. The AD brain networks seem to be more vulnerable to the removal of multiple types of couplings. A significant increase in gV in AD cases is observed in δ-γ, θ-β, α-β, β-γ, their inverse couplings and γ-γ. Interestingly, removal of WFC causes generally a larger increase in vulnerability compared to CFC, suggesting that while CFC plays a crucial role in the brain networks, WFC seems to be dominant in the brain networks. However, it is worth noting that the significant gV differences were observed almost exclusively in CFC (with the exception of γ-γ).
Finally, lV is used to assess the vulnerability of segregation of the network to the removal of a certain coupling type. Numerical results of comparing lV are reported in Table I and visualised in Fig. 8. γ-γ is the most robustly linked to segregation which fits well with the evidence of high-frequency oscillations being related to local processing [41]. Moreover, this coupling is significantly more vulnerable in HC, which is likely related to the decreased γ activity in AD [5]. Also, removal of α-α causes a significant increase of lV in HC and this function is likely disrupted in AD. On the other hand, β-γ, γ-α and γ-β removal cause a significant increase in lV in AD cases. This suggests, that despite the decreased activity of highfrequency components in AD networks, high-frequency CFC still enables the network segregation.

D. Conclusions and Future work
We have demonstrated that CBS and CS detect similar differences between AD and HC networks, but CBS has an advantage over CS by including cross-frequency and nonlinear interactions. We report several significant differences in CFC both globally and on a node-level, suggesting that including CFC in a graph-theoretic analysis of brain networks is crucial to obtain a more detailed insight into their structure and function. Furthermore, we show that multilayer network analysis provides a simple yet powerful framework for representing and analysing the role of CFC in brain networks. We present a novel approach to elucidate the roles of different frequency components of EEG signals using this framework.
CFC had been suggested to be related to modulatory activity, i.e. slow band modulating the activity of fast oscillations. However, it remains unclear why lowfrequency CFC would be increased in AD and requires further in-depth study.
Next, although (cross-)bispectrum was shown to be a powerful tool to detect various types of WFC and CFC, such as phase-phase or phase-amplitude, CBS seems to capture an unknown mixture of these types of couplings. Therefore, a combination of bispectrum with other types of CFC methods might be a plausible direction for future research.
Furthermore, by relying on traditional frequency bands to define the layers of the networks, our framework might miss some CFC occurring on finer scales, e.g. interaction within one band. However, considering the CFC within only a few bands allows us to construct multilayer networks with a relatively small number of layers. Thus, we argue that relying on the five bands is necessary to introduce the CFC into network analysis without increasing the complexity significantly.
Finally, the presented multilayer network analysis focused only on how dependent or vulnerable the networks are on different types of frequency coupling to enable integration and segregation properties. Although these two properties are hypothesised to be crucial in brain networks, their analysis might be insufficient to elucidate the function of the frequency couplings. Thus, we suggest to focus on other graph-theoretic measures beyond integration and segregation in future work. Moreover, it might be possible to use the multilayer networks features to classify the EEG networks automatically. We speculate that including the information about CFC might lead to machine learning models that perform significantly better than models that use only WFC-based features. However, we do not explore this route in this work. Instead, we focus on the novel multilayer network analysis to show the importance of considering CFC in addition to WFC when reconstructing EEG functional connectivity networks.

ACKNOWLEDGEMENT
The original data collection was funded by a grant from the Alzheimer's Research UK (ARUK-PPG20114B-25). The views expressed are those of the author(s) and not necessarily those of the NHS, the NIHR or the Department of Health.