ASSESSMENT OF CROSS-FREQUENCY PHASE-AMPLITUDE COUPLING IN NEURONAL OSCILLATIONS By Tamanna Tabassum Khan Munia A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Electrical Engineering — Doctor of Philosophy 2021 ABSTRACT ASSESSMENT OF CROSS-FREQUENCY PHASE-AMPLITUDE COUPLING IN NEURONAL OSCILLATIONS By Tamanna Tabassum Khan Munia Oscillatory activity in the brain has been associated with a wide variety of cognitive processes including decision making, feedback processing, and working memory control. The high temporal resolution provided by electroencephalography (EEG) enables the study of variation of oscillatory power and coupling across time. Various forms of neural synchrony across frequency bands have been suggested as the mechanism underlying neural binding. Recently, a considerable amount of work has focused on phase-amplitude coupling (PAC)– a form of cross-frequency coupling where the amplitude of a high-frequency signal is modulated by the phase of low-frequency oscillations. The existing methods for assessing PAC have certain limitations which can influence the fi- nal PAC estimates and the subsequent neuroscientific findings. These limitations include low- frequency resolution, narrowband assumption, and inherent requirement of bandpass filtering. These methods are also limited to quantifying univariate PAC and cannot capture inter-areal cross- frequency coupling between different brain regions. Given the availability of multi-channel record- ings, a multivariate analysis of phase-amplitude coupling is needed to accurately quantify the cou- pling across multiple frequencies and brain regions. Moreover, the existing PAC measures are usually stationary in nature, focusing on phase-amplitude modulations within a particular time window or over arbitrary sliding short time windows. Therefore, there is a need for computation- ally efficient measures that can quantify PAC with a high-frequency resolution, track the variation of PAC with time, both in bivariate and multivariate settings and provide a better insight into the spatially distributed dynamic brain networks across different frequency bands. In this thesis, we introduce a PAC computation technique that aims to overcome some of these drawbacks and extend it to multi-channel settings for quantifying dynamic cross-frequency cou- pling in the brain. The main contributions of the thesis are threefold. First, we present a novel time- frequency based PAC (t-f PAC) measure based on a high-resolution complex time-frequency dis- tribution, known as the Reduced Interference Distribution (RID)-Rihaczek. This t-f PAC measure overcomes the drawbacks associated with filtering by extracting instantaneous phase and ampli- tude components directly from the t-f distribution and thus provides high resolution PAC estimates. Following the introduction of a complex time-frequency-based high resolution PAC measure, we extend this measure to multi-channel settings to quantify the inter-areal PAC across multiple fre- quency bands and brain regions. We propose a tensor-based representation of multi-channel PAC based on Higher Order Robust PCA (HoRPCA). The proposed method can identify the signif- icantly coupled brain regions along with the frequency bands that are involved in the observed couplings while accurately discarding the non-significant or spurious couplings. Finally, we intro- duce a matching pursuit based dynamic PAC (MP-dPAC) measure that allows us to compute PAC from time and frequency localized atoms that best describe the signal and thus capture the tempo- ral variation of PAC using a data-driven approach. We evaluate the performance of the proposed methods on both synthesized and real EEG data collected during a cognitive control-related error processing study. Based on our results, we posit that the proposed multivariate and dynamic PAC measures provide a better insight into understanding the spatial, spectral, and temporal dynamics of cross-frequency phase-amplitude coupling in the brain. Copyright by TAMANNA TABASSUM KHAN MUNIA 2021 To my mom and dad. . . v TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix LIST OF ALGORITHMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Cross-frequency Coupling (CFC) for Neuronal Oscillation . . . . . . . . . . . . . 2 1.2 Methods for Quantifying Phase Amplitude Coupling . . . . . . . . . . . . . . . . 5 1.2.1 Hilbert transform based PAC measure . . . . . . . . . . . . . . . . . . . . 5 1.2.2 Wavelet transform based PAC measure . . . . . . . . . . . . . . . . . . . . 7 1.3 Phase Amplitude Coupling Modulation Indices . . . . . . . . . . . . . . . . . . . 8 1.3.1 Mean-Vector Length (MVL) Modulation Index . . . . . . . . . . . . . . . 8 1.3.2 Direct Mean-Vector Length (dMVL) Modulation Index . . . . . . . . . . . 8 1.3.3 Phase Locking Value (PLV) Modulation Index . . . . . . . . . . . . . . . 9 1.3.4 Kullback–Leibler Modulation Index (KL-MI) . . . . . . . . . . . . . . . . 10 1.4 Background on RID-Rihaczek Time-Frequency Distribution . . . . . . . . . . . . 10 1.5 Organization and Contributions of this Thesis . . . . . . . . . . . . . . . . . . . . 13 Chapter 2 Time-Frequency Based Phase-Amplitude Coupling Measure . . . . . . . 16 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2 Time-Frequency Phase Amplitude Coupling (t-f PAC) . . . . . . . . . . . . . . . . 18 2.2.1 Illustration of the Proposed time-frequency based PAC . . . . . . . . . . . 20 2.2.2 Significance Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.3.1 tf-MVL for Synthesized data . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.3.2 Comparison of tf-MVL Comodulogram with Conventional Hilbert-MVL and wavelet-MVL Comodulogram . . . . . . . . . . . . . . . . . . . . . . 27 2.3.3 Comparison of Phase Amplitude Coupling Indices (tf-MVL, tf-PLV and tf-MI) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3.4 Comparison of Accuracy of tf-MVL with Hilbert-MVL . . . . . . . . . . . 32 2.3.5 Comparison of tf-MVL with wavelet-MVL for Stationary/non-stationary Coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.4 Evaluation on EEG Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 Chapter 3 Multivariate Analysis of Phase-Amplitude Coupling using Tensor Ro- bust PCA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 vi 3.2.1 Tensor Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.3.1 Time-Frequency PAC (t-f PAC) . . . . . . . . . . . . . . . . . . . . . . . . 51 3.3.2 Statistical Significance Testing . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.3 PAC Tensor Construction Across Frequency Bands and Subjects . . . . . . 53 3.3.4 PARAFAC based Multivariate PAC . . . . . . . . . . . . . . . . . . . . . . 54 3.3.5 HoRPCA based Multivariate PAC . . . . . . . . . . . . . . . . . . . . . . 55 3.4 Description of Datasets and Experiments . . . . . . . . . . . . . . . . . . . . . . . 57 3.4.1 Synthesized Multi-channel Data . . . . . . . . . . . . . . . . . . . . . . . 57 3.4.2 EEG Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.4.3 Synthesized Data Experiments . . . . . . . . . . . . . . . . . . . . . . . . 60 3.4.4 EEG Data Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.5.1 Synthesized Data Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.5.2 EEG Data Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.5.3 Interrelationship between Multivariate PAC and hormone level . . . . . . . 75 3.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 Chapter 4 Assessment of Dynamic Phase Amplitude Coupling Using Matching Pursuit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.2.1 Phase Amplitude Coupling (PAC) . . . . . . . . . . . . . . . . . . . . . . 84 4.2.2 Sliding Window based Dynamic PAC . . . . . . . . . . . . . . . . . . . . 85 4.2.3 Matching Pursuit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4.2.4 Matching Pursuit based Dynamic PAC (MP-dPAC) . . . . . . . . . . . . . 87 4.2.5 Synthesized Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.2.6 EEG Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.3.1 Synthesized Data Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.3.2 EEG Data Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 Chapter 5 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . . 107 5.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 vii LIST OF TABLES Table 3.1: Comparison of PAC detection power obtained with HoRPCA, gedCFC, PARAFAC and Averaging based methods for different SNR values. The table depicts the Mean±standard deviation of evaluation matrices for 20 repetition of each ex- periment. Here, Prec.: Precision, Rec.:Recall, F-meas.: F-measure. . . . . . . . 68 Table 3.2: Comparison of detected channel pair accuracy for various SNRs obtained with HoRPCA, PARAFAC and Averaging based methods. The table depicts the Mean±standard deviation of evaluation matrices for 20 repetition of each ex- periment. Here, Prec.: Precision, Rec.:Recall, F-meas.: F-measure. . . . . . . . 69 Table 3.3: Comparison of detected channel pair accuracy for variability in channel loca- tions obtained with HoRPCA, PARAFAC and Averaging based methods. The table depicts the Mean±standard deviation of evaluation matrices for 20 repeti- tion of each experiment. Here, Prec.: Precision, Rec.:Recall, F-meas.: F-measure. 70 Table 3.4: Comparison of detected channel pair accuracy for various number of coupled channel pairs obtained with HoRPCA, PARAFAC and Averaging based meth- ods. The table depicts the Mean±standard deviation of evaluation matrices for 20 repetition of each experiment. Here, Prec.: Precision, Rec.:Recall, F-meas.: F-measure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 Table 3.5: Classification of EEG error and correct responses based multivariate PAC net- works computed using averaging, gedCFC, PARAFAC and HoRPCA based methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 Table 3.6: Comparison of theta-gamma interchannel t-f PAC across different estradiol lev- els. Bold values indicate significant values after Bonferroni correction. Here, * indicates significantly different for α = 0.05 and ** indicates significantly different for α = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Table 4.1: The number of MP atoms per trial, electrodes, clusters along with the signifi- cantly coupled atoms per clusters for error and correct response types averaged across all electrodes and subjects. . . . . . . . . . . . . . . . . . . . . . . . . . 99 Table 4.2: The MP-dPAC averaged across all channels and subjects for error and correct response types for different time windows. . . . . . . . . . . . . . . . . . . . . 102 viii LIST OF FIGURES Figure 1.1: Types of Cross-Frequency Coupling. . . . . . . . . . . . . . . . . . . . . . . . 4 Figure 1.2: Phase Amplitude Coupling between high frequency ( fa ) amplitude compo- nent, A fa (t) and low frequency ( f p ) phase component φ f p (t). . . . . . . . . . . 4 Figure 1.3: Illustration of the Hilbert transform based phase amplitude coupling (PAC) procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Figure 1.4: Illustration of various phase-amplitude coupling indices. . . . . . . . . . . . . 11 Figure 2.1: Illustration of the computation of phase-amplitude coupling on synthesized data. 20 Figure 2.2: Representation of the synthesized data generated for the analysis. . . . . . . . 26 Figure 2.3: Comparison of proposed tf-MVL method with conventional PAC measures through comodulograms and Shannon entropy measure. . . . . . . . . . . . . 30 Figure 2.4: Performance comparison of the Phase Amplitude Coupling Measures tf-MVL, tf-PLV and tf-MI as a function of varying coupling strength. . . . . . . . . . . 31 Figure 2.5: Performance comparison of tf-MVL measure with Hilbert Transform method for synthesized dataset 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 Figure 2.6: The robustness comparison of proposed tf-MVL measure with Hilbert-MVL method for synthesized dataset 2. . . . . . . . . . . . . . . . . . . . . . . . . 36 Figure 2.7: tf-MVL phase and amplitude component representation of synthesized data. . . 38 Figure 2.8: Comparison of tf-MVL and Wavelet-MVL based comodulograms for station- ary coupling. The signal was generated by coupling the amplitude of fa = 50Hz by the phase of a f p = 7Hz signal. . . . . . . . . . . . . . . . . . . . . . 39 Figure 2.9: The performance of tf-MVL and Wavelet-MVL for non-stationary coupling. The signal was generated by coupling the amplitude of fa = 50Hz by the phase of a linear chirp with f p = 7Hz. . . . . . . . . . . . . . . . . . . . . . . . . . 40 Figure 2.10: Topo plots of tf-MVL values for (a) ERN and (b) CRN response types. . . . . 42 Figure 2.11: p-value topo-plots for ERN and CRN after error and correct response types respectively. The channels marked in pink showed significant difference be- tween ERN and CRN (Wilcoxon Signed Rank Sum Test with p < 0.05). . . . . 42 ix Figure 2.12: Comodulograms show the PAC between low and high frequency bands for (a) CRN and (b) ERN of EEG data for FCz channel. . . . . . . . . . . . . . . . . 43 Figure 3.1: Illustration of the proposed HoRPCA based multivariate t-f PAC measure. . . . 57 Figure 3.2: The projections of the dipole locations to the scalp for the phase and ampli- tude components of the synthesized data for the performance evaluation ex- periments with varying signal parameters. . . . . . . . . . . . . . . . . . . . . 63 Figure 3.3: Low-rank and sparse networks extracted for the error and correct responses averaged across subjects. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 Figure 3.4: PAC network of p-values showing significant difference between error and correct responses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Figure 3.5: The detected channel combinations across all subjects for (a) error response (Xer ); and (b) correct response (Xcr ) for theta frequency band. . . . . . . . . . 73 Figure 3.6: Theta-gamma PAC networks detected for (a) Error response and (b) Correct response. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 Figure 3.7: Box plot of the theta-gamma interchannel t-f PAC between FCz-POz channel pair across different hormone levels (N = 43). . . . . . . . . . . . . . . . . . 76 Figure 3.8: Scatter-plot depicting the relationship between the computed theta-gamma t-f PAC and estradiol values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Figure 3.9: Pearson correlation analysis for theta-gamma t-f PAC and estradiol values. . . 78 Figure 4.1: Illustration of the proposed matching pursuit based dynamic PAC method. . . . 90 Figure 4.2: Clustering performance of the extracted MP atoms from synthesized data us- ing the (tstart ), (tstop ) and (tspan ) parameters. . . . . . . . . . . . . . . . . . . 94 Figure 4.3: Comparison of Matching Pursuit (MP) based dynamic PAC approach with sliding window (SW) based dynamic PAC approach for synthesized data. . . . 95 Figure 4.4: Performance comparison of dynamic PAC values for MP-dPAC and sliding window based dPAC for different coupling durations. . . . . . . . . . . . . . . 96 Figure 4.5: Performance comparison of dynamic PAC values for MP-dPAC and sliding window based PAC for different noise levels. . . . . . . . . . . . . . . . . . . 97 Figure 4.6: Coupling strength over time MP-dPAC map for error (first row) and correct (second row) responses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 x Figure 4.7: Dynamic PAC for EEG data during error and correct response for FCz channel in terms of comodulogram and coupling strength over time map. . . . . . . . . 101 Figure 4.8: Dynamic PAC for EEG data during error and correct response for F6 channel in terms of comodulogram and coupling strength over time map. . . . . . . . . 102 Figure 4.9: p-value topo-plots showing significant difference between error and correct responses for theta-gamma (first row) and alpha-gamma (second row) coupling. 103 xi LIST OF ALGORITHMS Algorithm 1: Sparse Representation-based Classification (SRC) of PAC Networks . . . . . . . . . 67 xii Chapter 1 Introduction The human brain has been modeled as a complex network with distributed topology. This dis- tributed topology results in parallel and specialized information processing. Information process- ing requires a neural mechanism to enable information integration and binding across specialized brain regions. Various forms of neural synchrony between oscillations across different frequency bands have been suggested as the major mechanism of neural integration. Previous studies based on electrophysiological measurement of neural activity suggest that different frequency bands are responsible for distinct computational roles [1] as oscillations are thought to create synchronization across specialized brain regions to corroborate cognitive processing [2, 3]. The power and/or the synchronization measured within and across different frequency bands have been related to various cognitive and neuronal functions [4, 5, 6]. For example, the gamma band neuronal activity in the human brain has been demonstrated to play an important role in visual perception [7, 8], whereas alpha band oscillations in the occipital region were interpreted to be an indicator of reduced visual attention [1]. Within band synchronization of excitability at different frequencies are considered crucial for cognitive functions and have been proposed to serve both for “bottom-up” and ”top- down” processing of information in the brain [9]. The within-band synchronization of neuronal oscillations is measured using phase synchrony which quantifies the phase difference between pop- ulations [10]. While phase synchronization within frequency bands can achieve spatial integration, it does not explain how interactions between phase-synchronized networks of different frequencies 1 can be achieved. Most recently, emerging evidence from various studies suggests that oscillations from different frequency bands are not isolated and independent; consequently, they can interact with each other in the form of modulation [11, 12]. This interaction of oscillations across different frequency bands is referred to as “cross-frequency coupling” [13, 14, 15, 16] which is described in detail in the next section. 1.1 Cross-frequency Coupling (CFC) for Neuronal Oscillation CFC has been proposed to be the “indicator of network coordination and functional integration” [17] and assumed to control the “integration of distributed information” [15]. It is also hypothe- sized as the ”mechanism for selectively routing information through neuronal networks” [18]. In essence, CFC is the phenomenon that is thought to facilitate communication between neuronal oscillations, as well as, enable parallel processing by separating the spatially distributed networks from one another. CFC measured between the time series of the same neuronal population is referred to as local CFC, whereas CFC measured between the time series of different neuronal populations is called inter-areal CFC. Multiple forms of cross-frequency coupling have been proposed including phase/amplitude [14, 19, 20, 21]; phase/phase [22, 9, 23, 24]; amplitude-to-amplitude [25, 26]; and phase-frequency [13, 16], as presented in Figure 1.1. Phase-to-phase coupling determines the statistical dependence between the instantaneous phase of signal X and instantaneous phase of signal Y1 as shown in Figure 1.1(a) (the phase-to-phase relation is 1:3 in this case) [27, 28, 25]. Phase-to-phase coupling is quantified using the cross-frequency phase synchrony (CFS) measures and it is hypothesized to control the ”binding of distributed neuronal activity” [9, 15]. This is also limited to n:m fre- quency ratios. Amplitude-to-amplitude coupling calculates the statistical dependence between the 2 instantaneous amplitude of frequency X and instantaneous amplitude of frequency Y2 of either the same or two distinct signals as shown in Figure 1.1 (b) [29, 15]. Amplitude-amplitude coupling has been reported in neuroimaging studies [25, 30, 31], but its functional significance is unclear [14] since such coupling is independent of spike time interactions [32]. Power to frequency cou- pling quantifies the changes in the frequency produced by variations in the amplitude envelope of the signal (measure the changes between the envelope of signal X and frequency modulations of signal Y4 as shown in Figure 1.1 (c)) [15]. The power/amplitude to frequency coupling is the least explored until recently and very few studies have reported it [33, 34, 35]. This type of coupling cannot be identified by other CFC measures, but it seems to play a critical role in the highly dy- namic systems, which can only be explained by ”frequency/phase entrainment” [33]. Furthermore, amplitude-frequency couplings can be recognized as so called autoresonance, when alterations in the driver frequency produces a corresponding alteration in the oscillation envelope, which results into prolonged phase locking of the driving and the oscillator frequency [33]. It is also well estab- lished in literature that frequency coupling as compared with amplitude coupling is less suscep- tible to interference and allows a higher dynamic range of the information processing. However, these aspects of information processing have not been well-explored in neuroscience until recently. Although there have been analyses using intracranial recordings [34, 35], power to frequency cou- pling is challenging to analyze in EEG and MEG data, where instantaneous frequency peaks are hard to determine [16]. In contrast to other forms of CFC measures, PAC is the most commonly used and best studied type of CFC measure. As shown in Figure 1.1 (d), PAC quantifies the statistical dependence between the instantaneous phase of frequency X and the instantaneous amplitude of frequency Y4 of either the same or two different signals. For a signal x(t), PAC between the high frequency ( fa ) and low frequency ( f p ) oscillations is defined as the modulation of the high frequency amplitude, 3 Figure 1.1: Types of Cross-Frequency Coupling. Figure adapted from J., Viktor, and V. Müller. ”Cross-frequency coupling in real and virtual brain networks.” Frontiers in computational neuro- science 7 (2013): 78. A fa (t), by the phase of the low frequency, φ f p (t) as shown in Figure 1.2. Figure 1.2: Phase Amplitude Coupling between high frequency ( fa ) amplitude component, A fa (t) and low frequency ( f p ) phase component φ f p (t). Figure adapted from R. T. Canolty and R. T. Knight, ’The functional role of cross-frequency coupling,’ Trends in Cognitive Sciences, 2010. PAC is thought to be responsible for integration across populations of neurons [36, 10]. It is also hypothesized to control the information transfer between distinct neuronal sites in mammalian 4 brains. Phase-amplitude coupling has been detected in different species including EEG, MEG and ECoG recordings of the human brain [19, 25, 36, 37], LFP of mice [38, 39], rats [40], sheep [41], monkeys [42] and within various brain regions such as the hippocampus [37, 38, 39, 42], the neocortex [12, 42, 43], and basal ganglia [19]. The neural information processing and cognitive functioning, particularly sensory signal detection [44], attention [20], visual perception [45, 46] and memory processing [21, 37] have also been shown to exhibit dynamic PAC. Disruption in PAC patterns has been linked to various neurological disorders, such as, autism spectrum disorder [47], schizophrenia [48], and Parkinson’s Disease [49]. More recently, PAC has been shown to be associated with cognitive control studies [50, 51]. Although a selection of methods has been developed in literature, no single method has been selected as the gold standard for assessing PAC. Different measures have their own advantages and shortcomings and may be used to serve certain purposes. Given the diverse role of PAC in neuronal functioning, there is a need to obtain unbiased, robust estimates of PAC. Therefore, in this thesis, we focus on the assessment of PAC in neuronal oscillations. 1.2 Methods for Quantifying Phase Amplitude Coupling Methods for quantifying phase amplitude coupling include Hilbert transform based and wavelet transform based measures described as follows. 1.2.1 Hilbert transform based PAC measure The most commonly implemented PAC measure is Hilbert transform based measure. In this ap- proach, PAC is estimated through the following steps [52, 53, 54]: 1) bandpass filtering the input data into bands of interest, e.g., theta and gamma; 2) applying Hilbert transform to extract ampli- 5 tude and phase time series from each frequency band of interest; and 3) quantifying the relationship between the phase and the amplitude time series. These steps are schematically illustrated in Figure 1.3. Figure 1.3: Illustration of the Hilbert transform based phase amplitude coupling (PAC) pro- cedure. Figure adapted from R. A. Seymour, G. Rippon and K. Kessler, ’The detection of phase amplitude coupling during sensory processing,’ Frontiers in neuroscience, 11, 487, 2017. However, Hilbert transform based amplitude and phase estimates suffer from some short- comings which can influence the final PAC estimates and the subsequent neuroscientific findings [53, 55, 56]. Most of these limitations arise from the first step of the processing pipeline, i.e. bandpass filtering. It is well-known that the use of Hilbert transform based amplitude and phase estimates rely on the assumption that the signals are narrowband, i.e., nearly sinusoidal [53, 57]. The narrowband assumption is questionable for high frequency oscillations. Therefore, using the Hilbert transform for extracting amplitude time series may result in non-meaningful amplitude es- timates for PAC computation [58, 59, 60]. Moreover, while using bandpass filtering for Hilbert transform based methods, it is advisable to use a wideband filter for high frequency oscillations and a narrowband filter for low frequency oscillations [61]. Therefore, in conventional PAC com- putation, a systematic bias arises due to the selection of the various bandpass filter parameters, 6 including the bandwidth, order of the filter, and transition band. Prior research indicates that the bandwidth for the high frequency component should be at least twice as high as the low frequency component in order to capture the amplitude modulation effect [58, 61]. Therefore, the bandwidth is proportional to the center frequency, and this choice can lead to a systematic bias. Moreover, Hyafil et al. [59] also showed that certain bandwidth choices might mistakenly result in the com- putation of phase-frequency coupling as PAC or produce erroneous amplitude-amplitude coupling. 1.2.2 Wavelet transform based PAC measure More recently, generalized Morse wavelet (GMW) has been proposed as an alternative to Hilbert transform for extracting the amplitude and phase components [62, 63] to address these limitations. GMWs are analytic wavelets, thus convolving the wavelet with the signal converts the signal to its analytic representation and is equivalent to bandpass filtering with a filter whose bandwidth scales with frequency. In this manner, the wavelet transform based PAC measure avoids the problem of designing the best filter as the filter bandwidth is automatically controlled by the wavelet scale. Estimating PAC from GMWs proceeds in the same way after the instantaneous amplitude and phase are extracted from the wavelet transform. Even though the wavelet transform addresses the problem of bandpass filtering, it suffers from a number of other problems including the choice of different design parameters in GMW [63]. Moreover, the lower and upper frequency and the scale parameters need to be chosen carefully to make sure that the peaks of the wavelet’s frequency spec- trum occur at the frequency of interest. Therefore, if the user does not choose the input frequency range or the number of scales appropriately, the ridges may not be as apparent. Consequently, PAC estimation may suffer from this. 7 1.3 Phase Amplitude Coupling Modulation Indices After extracting the phase and amplitude components using the Hilbert/wavelet based methods, the modulation between the low frequency ( f p ) phase component φ f p (t) and high frequency ( fa ) amplitude component A fa (t) is quantified using a PAC modulation index. PAC modulation in- dices available in literature includes the Mean Vector Length (MVL) modulation index, originally described in Canolty et al. [12]; the amplitude normalized direct Mean Vector Length (dMVL) modulation index described in Özkurt and Schnitzler [64]; the phase-locking value (PLV) modu- lation index described in [65]; and the Kullback–Leibler modulation index (KL-MI) described in Tort et al. [52]. 1.3.1 Mean-Vector Length (MVL) Modulation Index MVL index quantifies PAC by combining the phase component φ f p (t) and amplitude component A fa (t) to generate a complex-valued signal, in which each complex value corresponds to a vector in the polar plane. If the resultant distribution is non-uniform, this indicates a coupling between f p and fa , which can be quantified by taking the length of the average vector. MVL is calculated as: 1 T MV L( f p , fa ) = ∑ A fa (t)e jφ f p (t) , (1.1) T t=1 where T is the total number of time points. 1.3.2 Direct Mean-Vector Length (dMVL) Modulation Index MVL index has been reported to be partially dependent on the absolute value of the amplitude providing oscillation. As such any amplitude outliers can strongly influence the MVL index. To 8 overcome this limitation, Özkurt and Schnitzler [64] proposed an amplitude normalized MVL index as follows: T ∑t=1 A fa (t)e jφ f p (t) 1 dMV L( f p , fa ) = √ q , (1.2) T T ∑t=1 faA (t) 2 where T is the total number of time points. As dMVL is amplitude normalized, the value of dMVL ranges between 0 and 1. 1.3.3 Phase Locking Value (PLV) Modulation Index PLV index quantifies the modulation between the phase component φ f p (t) and amplitude com- ponent A fa (t) assuming that if PAC is present, the envelope of high frequency oscillation ( fa ) will oscillate at the frequency corresponding to low frequency oscillation ( f p ). The phase of high frequency envelope can be obtained from the angle of the analytic amplitude component. The modulation between the phase component φ f p (t) and amplitude component A fa (t) can be quan- tified by computing the circular variance of the consistency of the phase differences between the phase of the low frequency signal and the phase of the amplitude of the high frequency signal [66] as follows: 1 T j(φ f p (t)−φ fa (t)) PLV ( f p , fa ) = ∑e , (1.3) T t=1 where T is the total number of time points and φ fa (t) is the phase of high frequency envelope. 9 1.3.4 Kullback–Leibler Modulation Index (KL-MI) Kullback–Leibler modulation index (MI) quantifies PAC by calculating the deviation of the phase- amplitude distribution from the uniform distribution using Kullback-Leibler divergence [52]. First, a histogram is constructed from the phase angles in the range of -180 to 180 using an arbitrarily selected number of bins. The mean amplitude of the amplitude providing high frequency ( fa ) oscillation within each phase bin of the low frequency ( f p ) oscillation, is computed and normal- ized by the average value across all bins. The modulation index is quantified by comparing the phase-amplitude distribution (P) in contrast to the null hypothesis of a uniform phase-amplitude distribution (Q) as follows: D(P, Q) MI( f p , fa ) = . (1.4) log(Nbins ) The deviation between the two distributions is quantified by computing the Kullback-Leibler distance, related to Shanon entropy as follows: Nbins bin )  P(i  D(P, Q) = ∑ P(ibin )log . (1.5) ibin =1 Q(ibin ) 1.4 Background on RID-Rihaczek Time-Frequency Distribu- tion Rihaczek distribution is a complex time-frequency distribution [67] defined as: C(t, f ) = x(t)X ∗ ( f )e− j2π f t . (1.6) 10 Figure 1.4: Illustration of various phase-amplitude coupling indices. Figure adapted from M. J. Hulsemann, E. Naumann, and B. Rasch, “Quantification of phase-amplitude coupling in neuronal oscillations: Comparison of phase-locking value, mean vector length, modulation index, and generalized linear modeling cross-frequency coupling,” Frontiers in neuroscience, vol. 13, pp. 573, 2019. where x(t) is the signal and X( f ) is its Fourier transform. C(t, f ) measures the complex energy of a signal around time t and frequency f . The complex energy density function offers a deeper understanding of the characteristics of phase-modulated signals that is not obtainable with other time-frequency distributions. In contrast to wavelet transform, where the time-frequency resolution is determined by the wavelet function that expand the signal, for the Rihaczek distribution, the time-frequency resolution is quantified by the rate of change of the instantaneous frequency. As such, a better localization for phase- modulated signals are obtained. Similar to other members of Cohen’s class of distributions, the Rihaczek distribution is a bilinear time-frequency distribution and time and frequency shift covari- ant. This distribution satisfies the time and frequency marginals and preserves the signal energy with strong time and frequency supports [67]. With these properties, the Rihaczek distribution is identified as a complex time-frequency distribution that offers both time-varying energy spectrum 11 and the phase spectrum with good time-frequency localization [67]. One of the shortcomings of the Rihaczek distribution is the presence of cross-terms for multi- component signals. To overcome this limitation, the Reduced Interference Rihaczek (RID-Rihaczek) distribution is introduced. RID-Rihaczek distribution is an updated version of the Rihaczek distri- bution that employs the Choi-Williams kernel to filter out the cross-terms in the original distribu- tion. The distribution is defined as follows [68]: (θ τ)2 ZZ     θτ C(t, f ) = exp − exp j A(θ , τ)e− j(θt+2π f τ) dτdθ , (1.7) σ 2 2 where exp ( j θ2τ ) corresponds to the kernel function for the Rihaczek distribution [67], exp (− (θστ) ) corresponds to the Choi-Williams kernel, and A(θ , τ) is the ambiguity function of the given signal x(t) defined as: Z τ τ A(θ , τ) = x(u + )x∗ (u − )e jθ u du. (1.8) 2 2 This distribution still belongs to Cohen’s class of distributions and reflects both the time- varying energy and the phase of the signal [68]. The kernel functions can be thought of as two- dimensional lowpass filters that act on the ambiguity function that captures the time-varying au- tocorrelation of the signal. Through the kernel function, it is possible to reduce the effect of the cross-terms and localize the energy and phase estimates. RID-Rihaczek distribution has been previously employed for estimating phase and computing phase synchrony within a frequency band [68, 69]. It has been shown that RID-Rihaczek based phase synchrony measure has lower bias, better robustness to noise and higher time-frequency resolution compared to conventional methods including the Hilbert and wavelet transforms [68, 69]. 12 1.5 Organization and Contributions of this Thesis Although conventional methods for quantifying PAC have contributed greatly to the study of cross- frequency PAC, these measures present some drawbacks as mentioned in section 1.2. In this thesis, we present a novel PAC assessment technique that aims to overcome some of these drawbacks and extend the proposed method in some key ways to address the challenges with multichannel, non- stationary and directional couplings. In Chapter 2, we introduce a novel method for assessing PAC based on a high resolution com- plex time-frequency distribution. First, we introduce an instantaneous phase and amplitude es- timation method based on RID-Rihaczek distribution [68, 69]. The properties of this quadratic time-frequency distribution, such as time marginal and energy preservation, are used to estimate both the envelope of the high frequency oscillation and the phase of the low frequency oscillation. This approach replaces the Hilbert transform or the analytical wavelet transform for extracting the amplitude and phase. Unlike Hilbert transform based methods, the proposed method obtains the analytic signal without any bandpass filtering. Unlike the wavelet transform, this approach re- sults in uniformly high resolution across time and frequency and does not depend on the choice of different input parameters. We then compute MVL based on the extracted amplitude and phase time series to quantify PAC. It is important to note that even though the current chapter focuses on MVL to quantify the PAC, it is possible to combine the amplitude and phase estimates obtained from RID-Rihaczek distribution using other metrics such as MI or PLV. The proposed method is evaluated in terms of its resolution, accuracy of estimating the coupling strength and robustness against varying signal parameters for both simulated and real EEG signals. In Chapter 3, we extended our previously proposed time-frequency based PAC measure to quantify inter-areal PAC in multichannel recordings that can identify the low-high frequency oscil- 13 lation pairs with significant inter-areal coupling and determine their spatial locations. To accom- plish this, we employ a tensor representation of PAC across all channels, frequency bands, and subjects, to quantify the dynamics of multichannel CFC. Multi-way analysis of pairwise PAC in EEG takes full advantage of the multilinear structure of the data to provide a better understanding of the interactions across multiple modes. To that end, we propose a Higher order Robust Principal Component Analysis (HoRPCA) based multivariate framework to capture the spatial and spectral variation of the event-related PAC by employing the tensors constructed from PAC values across different channels, frequency bands, and subjects. To evaluate the efficacy of the proposed mea- sure, we assess its performance on simulated EEG data with varying signal parameters such as noise, subject variability, and volume conduction and compare it with existing multivariate PAC approaches. Finally, we apply the proposed method to EEG recordings collected during a cognitive control study to identify the spatial and spectral dynamics of pairwise PAC associated with differ- ent response types and to differentiate between the response types using multivariate analysis of PAC networks. Our results demonstrate that the proposed multivariate phase-amplitude coupling method can capture the spatial and spectral dynamics of PAC more accurately compared to ex- isting methods. Accordingly, we posit that the proposed higher order robust principal component analysis based approach filters out the background PAC activity and predominantly captures the event-related PAC dynamics to provide insight into the spatially distributed brain networks across different frequency bands. In Chapter 4, we focus on extending the proposed PAC measure to consider the dynamic na- ture of PAC. We introduce a PAC assessment technique that can quantify the temporal variation of PAC. Current methods for calculating PAC provide a numerical index representing an average value across a pre-determined time window. However, there is growing empirical evidence that PAC is dynamic, varying across time. Time-varying PAC measures available in the literature rely 14 on computing PAC over sliding short time windows. These approaches do not adapt to the sig- nal dynamics, and thus the arbitrary selection of the window length substantially hampers PAC estimation. To address these limitations, in this chapter, we introduce a data-driven dynamic PAC assessment approach to quantify time-varying PAC. The proposed approach relies on decompos- ing the signal using matching pursuit (MP) to extract time and frequency localized atoms that best describe the given signal. These atoms are then used to compute PAC across time and frequency. As the atoms are time and frequency localized, we only calculate PAC across time and frequency regions determined by the selected atoms rather than the whole time-frequency range. We evalu- ate the proposed method on both synthesized and actual EEG data. The results from synthesized data show that the proposed method detects the coupled frequencies and the time variation of the coupling correctly with high time and frequency resolution. The analysis of EEG data revealed theta-gamma and alpha-gamma PAC during response and post-response time intervals. As such, we hypothesize that the proposed MP based data-driven approach offers a more robust and possibly more sensitive method to quantify and track dynamic PAC effectively. Finally, in Chapter 5, we provide a summary of the dissertation and outline directions for future work. 15 Chapter 2 Time-Frequency Based Phase-Amplitude Coupling Measure Note: Most of the parts of this chapter have been published in the following article: 1. ”Munia, Tamanna TK, and Selin Aviyente. ”Time-frequency based phase-amplitude coupling measure for neuronal oscillations.” Scientific reports 9, no. 1 (2019): 1-15.” Note: Some parts of this chapter have been published in the following article: 2. ”Munia, Tamanna TK, and Selin Aviyente. ”Comparison of wavelet and rid-rihaczek based methods for phase-amplitude coupling.” IEEE Signal Processing Letters 26, no. 12 (2019): 1897- 1901.” 2.1 Introduction Neuronal oscillations across different frequencies play an important role in motor and cognitive functioning [70, 71]. Many ongoing hypotheses suggest that the coupling across frequencies, known as cross-frequency coupling (CFC), controls multi-scale information processing [72]. The most commonly studied type of CFC looks at the coupling between the amplitude of a high fre- quency oscillation and the phase of a low frequency oscillation and is known as Phase-Amplitude Coupling (PAC) [36, 10]. PAC between the amplitude of broadband gamma activity (30-100 Hz) 16 and the phase of various low frequency rhythms (typically 5-12 Hz) has been reported across sev- eral regions of the brain including the hippocampus, the basal ganglia, and the neocortex [71]. The frequency bands, the magnitude of coupling and the phase involved in PAC vary with time and with anatomical specificity during the execution of different cognitive and sensory tasks [14]. Since PAC plays an important role in quantifying nonlinear interactions between oscillatory processes, it is crucial to obtain an unbiased, robust estimate of PAC. Conventional PAC measures mostly rely on first bandpass filtering the data, followed by the Hilbert transform of filtered oscilla- tions to obtain phase and amplitude estimates of the oscillations [52]. However, these oscillations are not necessarily stationary or narrowband. Therefore, Hilbert transform based PAC computation may lead to erroneous assessment of PAC due to a systematic bias that can arise from the choice of the bandwidth, order of the filter, and transition band [73, 74, 75]. More recently, analytic wavelets such as Morlet and Generalized Morse Wavelet (GMW) have been employed to quan- tify PAC [62]. Even though the wavelet transform addresses the problem of bandpass filtering; it suffers from other problems including the dependence of the PAC estimates on the selected range of scales, the non-uniform time-frequency resolution and the different design parameters in GMW [62]. In this chapter, we introduce a novel approach for assessing PAC based on a high resolution complex time-frequency distribution. First, we introduce an instantaneous phase and amplitude estimation method based on RID-Rihaczek distribution [68, 69]. The properties of this quadratic time-frequency distribution, such as time marginal and energy preservation, are used to estimate both the envelope of the high frequency oscillation and the phase of the low frequency oscillation. This approach replaces the Hilbert transform and the analytical wavelet transform for extracting the amplitude and phase. Unlike Hilbert transform, the proposed method does not require bandpass filtering the signal to obtain the analytic signal. Unlike wavelet transform, this approach results 17 in uniformly high resolution across time and frequency and does not depend on the choice of different input parameters. We then compute MVL based on the extracted amplitude and phase time series to quantify PAC. Then we provide an analytical proof-of-concept for the RID-Rihaczek [76] time-frequency distribution based PAC measure. The proposed method is tested on simulated data and compared with Hilbert-based method in terms of resolution, accuracy of estimating the coupling strength and robustness against varying signal parameters. The method is also compared with the existing wavelet-based PAC measures for both stationary and non-stationary coupling mechanisms. Finally, the proposed method is applied to multi-channel EEG data collected from a cognitive control study to determine differences between response types and to identify brain regions and frequency bands that show increased PAC. 2.2 Time-Frequency Phase Amplitude Coupling (t-f PAC) For a signal x(t), with a high frequency component fa , and a low frequency component f p , PAC is defined as the modulation of the high frequency amplitude A fa (t), by the phase of the low frequency, φ f p (t). The first step is to extract the envelope of the high frequency amplitude signal and the phase of the low frequency signal. To compute PAC between low and high frequencies, first RID-Rihaczek distribution Cx (t, f ) is computed as in (1.7). To compute A fa (t), we propose to compute the frequency-constrained time marginal of Cx (t, f ) as follows: Zfa2 Axfa (t) = Cx (t, f )d f , (2.1) fa1 where fa1 and fa2 define the bandwidth around the high frequency of interest, fa . If the high frequency amplitude of the signal is coupled to the phase of a low frequency os- 18 cillation at f p , then this time marginal A fa (t) is expected to generate a peak at f p . This detected frequency is selected as the low frequency that is coupled with the amplitude of the high fre- quency component and the phase information at that frequency is extracted from the complex time-frequency distribution as follows:   C(t, f p ) φ f p (t) = arg . (2.2) |C(t, f p )| After detecting the amplitude and phase, we can use any existing PAC index such as mean vec- tor length (MVL) [12], Kullback-Leibler divergence modulation index (MI) [52] or phase locking value (PLV) [66]. MVL was suggested to be more sensitive to coupling strength and thus reported to be more suitable for high SNR data by Hulsemann et al. [54]. Therefore, in this chapter, the coupling between f p and fa is quantified by using MVL. This approach estimates PAC from a signal of length N, by mapping phase time series φ f p (t) given by (2.2) and amplitude time series A fa (t) given by (2.1) to a complex-valued vector at each time point, t [12]. To quantify the cou- pling between f p and fa , MVL measures the length of the average vector and computes PAC [12] as follows 1 N MV L( fa , f p ) = ∑ A fa (t)e jφ f p (t) . N t=1 (2.3) As MVL PAC index is used, for the remaining of this chapter, proposed t-f PAC method is referred as tf-MVL. The graphical representation of the proposed PAC computation approach is represented in Figure 2.1. 19 Figure 2.1: Illustration of the computation of phase-amplitude coupling on synthesized data. This synthesized data was generated by modulating the amplitude of an 80Hz high frequency signal with the phase of a 6Hz low frequency signal. (a) The synthesized signal; (b) Time-frequency component at high frequency; (c) Time marginal of the high frequency signal; (d) Time-frequency component at low frequency; (e) Illustration of detection of peak at coupled frequency from power spectral density of high frequency time marginal; (f) Phase component of the detected coupled low frequency; (g) t-f PAC coupling strength measurement using MVL method. 2.2.1 Illustration of the Proposed time-frequency based PAC For the illustration of the proposed measure, the sum of two sinusoids with known coupling param- eters was generated. These two sinusoids referred to as phase and amplitude signals are generated as follows: Phase signal : x f p (t) = K f p e j2π f pt , (2.4) 20 Amplitude signal : x fa (t) = A fa (t)e j2π fat , (2.5) where, A fa (t) = K fa (1 + m cos(2π f pt)), (2.6) m ∈ [0, 1], K f p and K fa are fixed scalars (for all simulations we set these to 1) that determine the amplitude of the phase frequency ( f p ) and amplitude frequency ( fa ), respectively. m determines the amount of modulation, thus the coupling strength. These two components are added to generate a synthesized signal: x(t) = x f p (t) + x fa (t). (2.7) A fa (t) and φ f p (t) extracted from the Rihaczek distribution are given below. These results can be easily extended to RID-Rihaczek distribution through a convolution with the kernel function. Given that x(t) = e j2π f pt + (1 + m cos(2π f pt))e j2π fat , (2.8) first X ∗ ( f ) is computed as follows: 1 1 X ∗ ( f ) = δ [( f − f p )] + mδ [ f − ( fa + f p )] + mδ [ f − ( fa − f p )] + δ [( f − fa )]. (2.9) 2 2 The Rihaczek distribution C(t, f ) is computed in (2.10) following (1.6) 21 1 1 C(t, f ) = δ [( f − f p )] + mδ [ f − ( fa + f p )]e− j2π fat + mδ [ f − ( fa − f p )] 2 2 1 e j2π(2 f p − fa )t + δ [( f − fa )]e j2π( f p − fa )t + δ [( f − f p )]e j2π( fa − f p )t + m 2 1 δ [ f − ( fa + f p )]e− j2π f pt + mδ [ f − ( fa − f p )]e j2π f pt + δ [( f − fa )] + m (2.10) 2 1 cos(2π f pt)δ [( f − f p )]e j2π( fa − f p )t + m2 cos(2π f pt)δ [ f − ( fa + f p )]e− j2π f pt 2 1 + m2 cos(2π f pt)δ [ f − ( fa − f p )]e j2π f pt + m cos(2π f pt)δ [( f − fa )]. 2 Extracting the time-marginal of C(t, f ) at f = fa provides the amplitude envelope as follows: Â fa (t) = (1 + m cos(2π f pt)). (2.11) The time-varying phase estimate based on the Rihaczek distribution can be derived as x(t)X ∗ ( f p )e− j2π f pt   φ̂ f p (t) = arg = φ f p (t). (2.12) |x(t)X ∗ ( f p )e− j2π f pt | Therefore, from (2.11) and (2.12), the amplitude estimate (Â fa (t)) and phase estimate (φ̂ f p (t)) obtained from Rihaczek distribution are equal to the actual amplitude (A fa (t)) and phase (φ f p (t)) around the frequency of interest. 2.2.2 Significance Testing The statistical significance of the tf-MVL measure was determined by generating surrogate datasets. Following the guidelines in Aru et al., we generated surrogate datasets using a block swapping ap- proach [55]. This method ensures that phase distortion is minimized, reducing the false-positive rate. A permuted signal was generated by splitting the envelope of the high frequency oscillation, 22 A fa (t) at a random time point and then swapping the two resulting time series. MVL values are then computed by associating the permuted A fa (t) with the original phase of the low frequency oscillation φ f p (t). This procedure of block swapping the amplitude time series and matching it with the original phase time series was performed 200 times to obtain the surrogate MVL values. This approach is known to produce relatively conservative evaluations of statistical significance of PAC measures [55]. This procedure preserves the global statistical characteristics of the data while disordering the correspondence between phase and amplitude time series and generating a sequence of MVLs observed under the null hypothesis that phase-amplitude coupling is due to chance. The observed tf-MVL value was deemed to be significant only when it is larger than 95% of the surrogate tf-MVL values. As the proposed time-frequency measure uses a quadratic distribution, the phase and amplitude estimates may be non-local causing a bias in the MVL estimates. However, these non-localities are minimized through the Choi-Williams kernel used in the definition of RID-Rihaczek distribution in (1.7). Moreover, Cohen et al. [57] assert that the bias in MVL estimates due to non-stationarities in the signals can be mitigated by applying shuffling based permutation or debiasing terms [77, 78]. 2.3 Results To investigate the validity of the proposed PAC approach, experiments were first conducted on two sets of synthesized data, and then on human EEG data taken from a previously published study [79]. 23 2.3.1 tf-MVL for Synthesized data Synthesized Dataset 1 In the first simulation, the sum of two sinusoids with known coupling parameters was generated to assess the accuracy of the proposed method. These two sinusoids referred to as phase and amplitude signals, are generated following Tort et al. as follows [52, 80]: Phase signal : x f p (t) = K f p sin(2π f pt), (2.13) Amplitude signal : x fa (t) = A fa (t) sin(2π fat), (2.14) (1−χ) sin(2π f p t−φc )+χ+1 where, A fa (t) = K fa 2 , χ ∈ [0, 1], K f p and K fa are constants that determine the amplitude of the phase frequency ( f p ) and amplitude frequency ( fa ), respectively. φc is the phase of the low frequency phase providing signal at the time point where the magnitude of the amplitude signal bursts is maximum. χ is the fraction of amplitude signal that is not modulated by the phase signal, thus the coupling strength is given by (1 − χ). To generate the synthesized dataset 1, these two components were added: x(t) = x fa (t) + x f p (t) + ε(t), (2.15) where ε(t) is additive noise and was generated with a combination of random samples created through power law and white Gaussian noise. The power law samples simulate the background brain activity, and the white Gaussian samples simulate measurement noise space [80]. The strength of the white Gaussian noise parameters was set to half of the strength of the power law samples. The synthesized dataset was generated with the low frequency signal at 5 Hz, high frequency signal in the range of 50-70 Hz, sampling frequency of 1000 Hz, SNR levels ranging from -5 dB to 24 10dB, coupling intensity in the range of 0.1 to 1.0, and data length in the range of 0.4 to 6 seconds. Finally, to test the effect of coupled high frequency, a 10 second signal was generated with SNR = 3 dB, sampling frequency = 1000 Hz, coupling intensity = 0.9, low frequency = 5 Hz, and changing the coupled high frequency from 10 Hz to 70 Hz with an increment of 5 Hz. Synthesized Dataset 2 The second synthesized dataset was generated using the fact that the spectrum of EEG data fol- lows the power law, i.e., higher the frequency, weaker the amplitude (P( f ) ∼ ( f1β )). The strength of the amplitude decrease is defined by the parameter β with β = 0, 1, and 2 indicating white noise, pink noise, and Brownian (red) noise, respectively [54]. The spectrum of EEG data is reported to be related to Brownian noise [81, 82], thus, the synthesized data was generated from Brownian noise. First, 10 seconds of Brownian noise data was generated at a sampling frequency of 1000 Hz following the method developed by Zhivomirov [83]. Next, the signal was bandpass filtered at low phase providing frequency with a bandwidth of 2 Hz for generating the phase signal. The same simulated signal was then bandpass filtered at high amplitude providing frequency to create the amplitude signal. The low frequency phase signal bandwidth was set to 8-10 Hz, and the high frequency amplitude signal bandwidth was set to 50-70 Hz. The coupling between phase and am- plitude signals was generated using the procedure described by Kramer and Eden [84]. The time locations of relative maxima and minima of the phase signal were detected. At each maxima, a DC shifted Hanning window with a duration of 42 ms, i.e., the amplitude of the Hanning window is shifted by 1, was multiplied with the amplitude time series. The monophasic coupling was gen- erated by multiplying the Hanning window with the amplitude time series centered at the relative maxima of the phase time series. The intensity of the phase-amplitude coupling is controlled by multiplying the Hanning window with a constant I, where I = 1 indicates full coupling and I = 0 indicates no phase-amplitude coupling. An additional frequency modulated noise was generated 25 by generating Brownian noise of the same length, bandpass filtering the noise signal at similar low and high frequency bands and adding them to the phase and amplitude signals, respectively. Synthesized data were generated with low frequency signal at 10 Hz, high frequency signal in the range of 50-70 Hz, sampling frequency of 1000 Hz, SNR levels ranging from -5 dB to 17dB, coupling intensity in the range of 0.1 to 1.0, and data length in the range of 0.4 to 6 seconds. To test the effect of coupled high frequency, a 10 second signal was generated with SNR = 6 dB, sampling frequency = 1000 Hz, coupling intensity = 0.9, low frequency = 10 Hz, and changing the coupled high frequency from 15 Hz to 70 Hz with an increment of 5 Hz. Finally for visualizing the broadband vs narrow-band coupling effect, 10 second signals (SNR = 6 dB, sampling frequency = 1000 Hz, coupling intensity = 0.85) were generated with the f p frequency fixed at 5 Hz and varying the fa frequency range as [65-75], [60-80], [55-85], [50-90], [45-95] and [40-100] to generate 10 Hz to 60 Hz coupling bandwidths. Examples from the two synthesized data sets are shown in Figure 2.2. Figure 2.2: Representation of the synthesized data generated for the analysis. (a) Simulated signal generated through the method described in synthesized dataset 1 subsection; (b) Simulated signal generated through the method described in synthesized dataset 2 subsection. Both signals were generated for 2 seconds with a sampling frequency = 1000Hz, f p = 5Hz, fa = 70Hz, SNR = 3 dB and coupling intensity = 0.7. 26 2.3.2 Comparison of tf-MVL Comodulogram with Conventional Hilbert- MVL and wavelet-MVL Comodulogram We compare the performance of three different approaches, i.e. Hilbert transform, wavelet trans- form and the proposed time-frequency approach, for computing PAC as reflected by the MVL comodulograms. The comparison of these three methods was performed by evaluating the co- modulograms both qualitatively and quantitatively. Comodulogram is a 2-D map that indicates the strength of coupling, as measured by (2.3), between different oscillation frequencies, with the f p values along the x-axis, and fa values along the y-axis [52]. As the coupling strength is plotted as a function of low and high frequency, comodulogram maps the coupling strength for ( fa and f p ) and exhibits highest value for the combination with highest PAC. In all these cases, the comodu- lograms were constructed by considering a frequency step size of 1Hz for fa and 1Hz for f p . For the comparison, a synthesized signal was generated following the synthesized dataset 1 method described above with high frequency fa = 60 Hz, low frequency f p = 10 Hz, coupling coefficient = 0.8, and SNR = 6 dB. The proposed method is first compared with the conventional Hilbert transform based method [54, 77, 85, 86]. The filter designed for the Hilbert transform based approach was a fourth order bandpass Butterworth filter. The bandwidth of the bandpass filter is selected to be variable as discussed in [55]. The filters that extract the envelope of the high frequency oscillation, fa , should be wide enough to accomodate the center frequency ± the modulating f p . If this condition is not met, the PAC cannot be identified without fulfilling this condition even if present [58]. As a result, a variable bandwidth defined as ±0.4 times the amplitude frequency (e.g. for fa = 40Hz, the chosen bandwidth was 40 ± 0.4 × fa = [24, 56]) was chosen, which has been shown to improve the detectability of PAC in literature for Hilbert transform [61, 86]. The bandwidth for f p was 27 kept narrow ( f p ± 1Hz) to extract the sinusoidal waveforms [86]. After filtering, Hilbert transform was performed to extract the instantaneous phase and amplitude components and finally PAC was calculated by computing the MVL for the Hilbert method (Hilbert-MVL). Generalised Morse wavelet (GMW) based PAC [62] was also computed for comparisons. This method generates complex valued time series through analytic wavelet transform from where phase and amplitude components can be extracted using this complex valued signal. As described in [62], the PAC measure using generalised Morse wavelet largely depends on the value of the selected parameters such as β and γ. In this chapter, we set the values of these parameters as β = 6 and γ = 3 as suggested in Nakhnikian et al. [62]. After extracting phase and amplitude components, we compute the MVL for generalized Morse wavelet (Wavelet-MVL) similar to tf-MVL and Hilbert- MVL. The resulting comodulograms are illustrated in Figure 2.3 (a-c). As it can be seen from this figure, the tf-MVL method provides higher resolution estimates of the two frequencies that are coupled with each other compared to Hilbert transform and generalized Morse wavelet trans- form based methods. It can also be seen that GMW based method has better localization of the frequencies compared to the Hilbert transform. After computing the comodulograms, the three methods were also compared quantitatively using Shannon entropy to quantify the concentration of the comodulograms. The Shannon entropy for a comodulogram is defined in the same way that entropy is quantified for time-frequency distributions [87], i.e. H = − ∑ f1 , f2 p f1 , f2 log2 p f1 , f2 , MV L( f1 , f2 ) where p f1 , f2 = ∑ f1 ∑ f2 MV L( f1 , f2 ) . The comodulogram surface is first normalized to obtain a two- dimensional distribution similar to a probability density function and then the definition of Shannon entropy is applied to this surface. In information theory, a maximum value of entropy is obtained for a uniform distribution, whereas the minimum value is achieved for a Dirac delta function, i.e. when there is no uncertainty about the value of the random variable. In a similar fashion, for a 28 comodulogram a high entropy value corresponds to a surface that is equally distributed across dif- ferent frequency values while a low entropy value corresponds to a surface that is well-localized. Therefore, in this chapter, entropy is used as a quantitative measure of concentration of the comod- ulogram surface similar to the way it has been applied to time-frequency distributions in prior work [87, 88]. For different SNR values, the Shannon entropy of the comodulograms is given in Figure 2.3 (d). The lower Shannon entropy value for tf-MVL indicates that the MVL modulation index is more concentrated, i.e., has a higher resolution, and thus leads to more accurate PAC detection. Moreover, the entropy for the tf-MVL method stays constant across different SNRs indicating the robustness of the method to noise. In comparison, wavelet and Hilbert transform based methods have higher entropy, with Hilbert transform yielding the highest entropy, i.e. the worst localization. Moreover, both the wavelet and Hilbert transform based methods show a fluctuation of the entropy values across SNRs, indicating less robustness compared to time-frequency based method. As the wavelet based method is fairly new and depends largely on the selection of the wavelet design parameters and the range of input frequencies, further research is required for the correct quantification of PAC using this approach. On the other hand, Hilbert transform based method is well-established. For these reasons, in the following sections, the proposed tf-MVL measure is compared to the conventional Hilbert-MVL approach in terms of the resolution of the comodulo- gram, accuracy of the estimated PAC value and robustness to varying signal parameters for both simulated and real data. 29 Figure 2.3: Comparison of proposed tf-MVL method with conventional PAC measures through comodulograms and Shannon entropy measure. (a) Hilbert transform based PAC (Hilbert-MVL) measure; (b) Generalized Morse Transform based PAC (Wavelet-MVL) measure and (c) Proposed time-frequency based PAC (tf-MVL) measure; (d) Comparison of Entropy Mea- sure for Hilbert-MVL, Wavelet-MVL and tf-MVL comodulograms for varying SNR levels. A low entropy value corresponds to a comodulogram that is well-localized, while a high entropy value corresponds to a comodulogram that is distributed across a wide range of frequencies. The en- tropy for the tf-MVL method stays constant across different SNRs indicating the robustness of the method to noise. 2.3.3 Comparison of Phase Amplitude Coupling Indices (tf-MVL, tf-PLV and tf-MI) To quantify the PAC, our proposed time-frequency based method can be used with any of the existing PAC indices such as MVL, PLV or MI, where the amplitude and phase are extracted using 30 the proposed time-frequency method. To compare the performance of different PAC measures, tf-MVL, tf-PLV and tf-MI were computed as a function of varying coupling strength. The model described in synthesized dataset 1 was used for this comparison. Data was generated with f p = 5 Hz, fa = 70 Hz, SNR = 6 dB with a sampling frequency of 1000 Hz. Five different coupling strengths ranging from 0.1 to 1 with a stepsize of 0.2 were considered. In all of these cases, RID-Rihaczek has been used to extract the amplitude and the phase of the high and low frequency components, respectively, but different metrics have been employed to quantify the PAC. As shown in Figure 2.4, the PAC values of all methods increase with increasing coupling strength. Figure 2.4: Performance comparison of the Phase Amplitude Coupling Measures tf-MVL, tf- PLV and tf-MI as a function of varying coupling strength. PAC values of all methods increased with increasing coupling strength and MVL differentiates best between different levels of coupling strengths. 31 In accordance with existing literature [54], MVL is found to be more sensitive to the varying coupling strength, followed by MI and PLV. Despite small quantitative differences in these three measures, the same trend is observed for all three methods in the simulation. For this reason, all of the following analysis was performed using tf-MVL method. 2.3.4 Comparison of Accuracy of tf-MVL with Hilbert-MVL Synthesized dataset 1 was used for this comparison as the ground truth for coupling strength is known. The data was generated with f p = 5 Hz and fa = 70 Hz with a sampling frequency of 1000 Hz. The range of interest for the high frequency was set to [40, 90] Hz whereas the low frequency range was chosen to be [2, 12] Hz. The two methods were compared by calculating the relative error rate of the estimated coupling strength using the following equation: PACdetected − PACactual RelativeError(RE) = × 100%. (2.16) PACactual Figure 2.5 shows the comparison of the proposed method with Hilbert transform in terms of quan- tification of coupling coefficient under four different conditions: varying SNR from -5 to 10 dB, varying coupling strength of 0.1 to 1, varying data length of 0.4 to 6 seconds and varying the high frequency of 10 Hz to 70 Hz. For varying SNR levels, it can be seen in Figure 2.5 (a) that the proposed method is more robust against noise with smaller errors in estimating the coupling parameter (MeanError < 10%). Even at high noise levels, e.g., SNR = −5 dB, the error rate is as low as 28% for tf-MVL compared to 63% relative error exhibited by Hilbert transform based method. For varying degrees of coupling strength, ten different coupling strengths ranging from 0.1 to 1 were considered. The relative errors in coupling strength estimation are depicted in Figure 2.5 32 (b). The proposed method is better at quantifying the true coupling strength compared to Hilbert transform based PAC. As tf-MVL does not depend on bandpass filtering, the performance of tf- MVL (Mean Error: 2.62 ± 3.27%) is much better at detecting the coupling strength compared to Hilbert Transform (Mean Error: 21.94 ± 11.36%). The sensitivity of PAC to data length is evaluated in Figure 2.5(c). Different signal lengths from 0.4 to 6 seconds were considered. For both methods, the accuracy of detecting the coupling strength improved with increasing signal length, although at different rates. Overall, for signal durations longer than one second, tf-MVL produced the best coupling strength estimates. Finally, the two methods were compared with respect to the difference between the two fre- quencies by changing the high frequency amplitude providing signal. The fa value was varied from 10 Hz to 70 Hz with an increment of 5 Hz by keeping the f p value constant at 5 Hz. As shown in Figure 2.5(d), the change of frequency has no effect on the proposed t-f based PAC measure result- ing in consistently low error rates, whereas Hilbert transform performed better as the difference between the two frequencies got larger. tf-MVL and Hilbert-MVL were compared for synthesized dataset 2 for different signal param- eters including varying SNR levels, coupling strength, data length, high frequency component and narrowband vs. broadband coupling. In this case, as there is no ground truth for the coupling strength, instead of the relative error, the two methods were compared using the MVL values. Moreover, the PAC values computed from the surrogate data for both the tf-MVL and Hilbert-MVL methods were compared with the respective observed MVL values to determine the significance of the results. As shown in Figure 2.6 (a), the tf-MVL method was more robust than the Hilbert-MVL for different noise levels. The tf-MVL method reached the actual PAC value at 1dB whereas the Hilbert-MVL method showed a high MVL value only at high SNR levels. Both tf-MVL and 33 Figure 2.5: Performance comparison of tf-MVL measure with Hilbert Transform method for synthesized dataset 1. (a) Relative error for various SNR levels; (b) Relative error for various coupling strengths; (c) Relative error for different data lengths; (d) Relative error for coupled high frequency. Hilbert-mvl showed an increase in MVL with an increase in coupling intensity. As shown in Figure 2.6 (b), tf-MVL method was more sensitive to the coupling strength compared to Hilbert- MVL. The effect of data length is shown in Figure 2.6 (c). The tf-MVL is more robust and reaches the level of acceptable PAC at 1.2 seconds compared to Hilbert-MVL which provides acceptable PAC only at 3.2 seconds, thus indicating that tf-MVL method is still effective when data length is very small. As shown in Figure 2.6 (d), the high frequency that the low frequency is coupled to has no effect on the tf-MVL method. As no filtering is required for tf-MVL, this method can detect 34 PAC even when the frequency difference between high and low frequency is as low as 5 Hz. On the other hand, due to the required filtering in the implementation of the Hilbert-MVL method, this method can only start to detect the true coupling when the high frequency is at 40 Hz. In Figure 2.6 (e), we evaluated the effect of narrowband vs. broadband coupling. For this comparison, the f p frequency was fixed at 5 Hz whereas the fa frequency range was varied as [65-75], [60-80], [55-85], [50-90], [45-95] and [40-100] to generate coupling bandwidth ranging from 10 Hz to 60 Hz. The SNR was fixed at 6dB with coupling strength fixed at 0.9. The variation of tf-MVL is low compared to Hilbert-MVL with varying coupling bandwidth. As expected, an increase in Hilbert-MVL was found with an increase in coupling bandwidth. The threshold values obtained through surrogate data testing (MVL threshold for 95% confi- dence interval) are also shown in Figure 2.6. As reflected in Figure 2.6 (a)-(e), for both tf-MVL and Hilbert-MVL cases the observed MVL was significantly higher (one tail Wilcoxon Signed Rank Sum Test, p < 0.05) than the threshold, thus rejecting the null hypothesis. 2.3.5 Comparison of tf-MVL with wavelet-MVL for Stationary/non-stationary Coupling In the first simulation, simulated data was generated following (2.8) by modulating the amplitude of high frequency signal at fa = 50Hz with the phase of a low frequency signal at f p = 7Hz with coupling strength m = 0.95. After generating the signal, the PAC value was computed first using tf-MVL and then GMW based PAC estimation method. Figure 2.7(a) shows the generated synthe- sized signal x(t). Figure 2.7(b) and (c) show the RID-Rihaczek time-frequency (t-f) distribution of the low and high frequency components of the signal. Figure 2.7(d) and (e) show the wavelet distribution of the low and high frequency components, respectively. The generated and estimated 35 Figure 2.6: The robustness comparison of proposed tf-MVL measure with Hilbert-MVL method for synthesized dataset 2. The broken lines indicate the critical threshold for PAC from 200 surrogate time series. (a) MVL for various SNR levels; (b) MVL for various coupling strengths; (c) MVL for various data lengths; (d) MVL for different coupled high frequency com- ponents; (e) MVL for various bandwidths of coupled high frequency. Synthesized signals were generated with the low-frequency phase component at 5 Hz and a varying frequency range for the high-frequency component as [65-75], [60-80], [55-85], [50-90], [45-95] and [40-100] to create PAC with 10 Hz to 60 Hz amplitude bandwidth. 36 magnitude and phase through both methods are shown in Figure 2.7(f) and in Figure 2.7(g), re- spectively. Comparing both methods, we can see that the proposed RID-Rihaczek method can correctly extract the phase and amplitude components whereas the wavelet method extracts a low frequency envelope estimate that does not capture the modulation effect while detecting the correct phase signal. The comodulogram for both tf-MVL and wavelet-MVL methods are shown in Figure 2.8(a) and (b) respectively to illustrate the strength of coupling between different oscillation frequencies. As exhibited in the figure, tf-MVL comodulogram has higher frequency resolution compared to wavelet-MVL comodulogram. For the next experiment, we generated non-stationary coupling by replacing the cosine function in Eq. (2.6) with a chirp as follows: A fa (t) = K fa (1 + m cos(2π f pt α )), (2.17) where α = 2. This non-stationary modulation is generated in order to explore the robustness of our method when non-stationary PAC is present in the system. The chirp signal was generated with a starting low frequency of f p = 7Hz and high frequency amplitude component at fa = 50Hz. Representative results for a linear chirp comparing both methods are shown in Figure 2.9. From the low frequency energy of the time-frequency distribution shown in Figure 2.9(b), we can see the modulating frequency ranges from 2-12 Hz for RID-Rihaczek and 2-20 Hz for wavelet. Therefore, the modulating phase component was extracted for the bandwidth 2-12 Hz for tf-MVL and 2-20 Hz for wavelet-MVL. Comparing both results in Figure 2.9(f) and Figure 2.9(g), we found that our method is more accurate in detecting the coupling between f p = 7Hz and fa = 50Hz. 37 Figure 2.7: tf-MVL phase and amplitude component representation of synthesized data. (a) Synthesized signal ( f p = 7 Hz and fa = 50 Hz; (b) Magnitude distribution at low frequency for tf-MVL method; (c) Magnitude distribution at high frequency for tf-MVL method; (d) Magnitude distribution at low frequency for wavelet-MVL method; (e) Magnitude distribution at high fre- quency for wavelet-MVL method; (f) Estimated envelope of the high frequency component of the synthesized signal; (g) Estimated phase of the synthesized signal. 38 Figure 2.8: Comparison of tf-MVL and Wavelet-MVL based comodulograms for stationary coupling. The signal was generated by coupling the amplitude of fa = 50Hz by the phase of a f p = 7Hz signal. The comodulogram shows the PAC between low frequency ( f p = 7Hz) and high frequency ( fa = 50Hz) components of the synthesized data for (a) tf-MVL based method and (b) for wavelet-MVL based method. 2.4 Evaluation on EEG Data The proposed PAC measure was evaluated on an EEG dataset from a previously published cogni- tive control-related error processing study [79]. A letter version of the speeded-reaction Flanker task [89] was designed following the experimental protocol approved by the Institutional Review Board (IRB) of the Michigan State University. The data collection was performed in accordance with the guidelines and regulation established by this protocol. Written and informed consent was collected from each participant before data collection. In this chapter, data from 19 participants were considered. A string of five letters, which could be congruent (e.g., SSSSS) or incongru- ent (e.g., SSTSS), was visually presented to each participant for each trial. The participants have to choose the center letter with a standard mouse with respect to the Flanker letters. The trials began with a 35ms of flanking stimuli (e.g., SS SS). The target stimuli were then embedded in the center of the flankers (e.g., SSSSS/SSTSS) and remained for 100 ms (total presentation time is 135ms) followed by an inter-trial interval of varying duration ranging from 1200 to 1700 ms. 39 Figure 2.9: The performance of tf-MVL and Wavelet-MVL for non-stationary coupling. The signal was generated by coupling the amplitude of fa = 50Hz by the phase of a linear chirp with f p = 7Hz. (a) Representative 1 second epoch of the signal; (b) Low frequency energy distri- bution for RID-Rihaczek; (c) Low frequency energy distribution for wavelet; (d) High frequency energy distribution for RID-Rihaczek; (e) High frequency energy distribution for wavelet; (f) Co- modulogram generated using tf-MVL; (g) Comodulogram generated using wavelet-MVL. The trials were conducted to study the Error-Related Negativity (ERN) after an error response and the Correct-Related Negativity (CRN) after a correct response. Total number of trials was 480 in which the total number of error trials in different participants varied from 20 to 61. To keep 40 the comparison between ERN/CRN consistent, the same number of correct responses was chosen randomly from the correct trials. Continuous EEG responses were recorded using the ActiveTwo system (BioSemi, Amsterdam, The Netherlands) by placing the 64 Ag–AgCl electrodes in accor- dance with the international 10/20 system. All EEG signals were sampled at 512 Hz. The trials with artifacts were removed, and the volume conduction was minimized using the Current Source Density (CSD) Toolbox [90]. The artifact removed data averaged over all trials were considered for the analysis in this chapter. The time-frequency based PAC measure was applied to EEG data collected during a cognitive control study where ERN and CRN were compared after error and correct responses, respectively. As previous studies indicate increased synchronization associated with the ERN in the time win- dow 25-75 ms [69, 79], all analysis was performed for the 25-75 ms time window of the EEG data. For each electrode, PAC was computed between a band of low frequency oscillations ( f p ) in the [2, 13] Hz range, and high frequency oscillations ( fa ) in the [34, 100] Hz range. The average cou- pling strength was calculated for 66 frequencies distributed linearly in the fa range of interest with an increment of 1 Hz. The topoplots showing the computed tf-MVL value for all 64 electrodes for both ERN and CRN are given in Figure 2.10. From this figure, it can be seen that the PAC averaged across subjects is higher for ERN compared to CRN and that the high PAC values for ERN are concentrated around FCz and central-parietal regions. A Wilcoxon Signed Rank Sum Test is conducted to determine the statistical significance of the difference between ERN and CRN PAC values across all electrodes. tf-MVL value is significantly higher (p = 0.00381) for ERN compared to CRN. A p-value topoplot that shows the channels with significant difference between ERN and CRN is given in Figure 2.11. The most significant difference was found for FCz with p = 0.000367 after Bonferroni correction. 41 Figure 2.10: Topo plots of tf-MVL values for (a) ERN and (b) CRN response types. Figure 2.11: p-value topo-plots for ERN and CRN after error and correct response types respectively. The channels marked in pink showed significant difference between ERN and CRN (Wilcoxon Signed Rank Sum Test with p < 0.05). 42 As FCz provided the highest tf-MVL value and the highest significant difference between the two response types, the coupling strength at FCz is depicted as a comodulogram in terms of fa vs. f p , for both error and correct responses in Figure 2.12. The overall strongest PAC index was found between the phase of alpha (9-12 Hz) and the amplitude of slow gamma (40-60 Hz) oscillations (Figure 2.12). Though both ERN and CRN responses exhibited high PAC for these ranges of frequencies, the coupling coefficient for ERN (MV L = 0.22424±0.1039) was significantly greater compared to the coupling coefficient for CRN (MV L = 0.0.06758 ± 0.03579) with a p = 0.0031 determined by a Wilcoxon Signed Rank Sum Test. Figure 2.12: Comodulograms show the PAC between low and high frequency bands for (a) CRN and (b) ERN of EEG data for FCz channel. 2.5 Discussion Human brain is known to route information through multiple networks operating in parallel. Os- cillatory coupling across frequency bands provides the temporal and spatial dynamics necessary to enable these networks. One prominent type of cross-frequency coupling known as phase-amplitude coupling may reflect integration processes across populations of neurons. PAC has been reported to control the long-distance communication considering that the slow oscillations of the neuronal sig- 43 nal can propagate at larger scales compared to the fast ones [13, 11, 16] and thus was proposed as ”canonical mechanism for neural syntax” [53]. Although there is a growing literature on measures to quantify PAC from neurological data, the best approach to detect and quantify the phenomena is still difficult to settle on. In this chapter, we proposed a time-frequency based PAC as a novel method for estimating cross-frequency phase-amplitude coupling and illustrate that the proposed method offers higher accuracy and robustness compared to existing methods through simulations. We first tested the model on two different sets of simulated signals, allowing us to generate a fully controlled com- parison of models and parameters and then applied the proposed tf-MVL model to estimate the PAC from neurophysiological signals. The results from synthesized data showed that the proposed time-frequency based PAC method provides several advantages over existing Hilbert and Wavelet- MVL methods. First, the proposed method estimates both the amplitude of the high frequency signal component and the phase of the low frequency component directly from a high resolu- tion complex time-frequency distribution. Using properties of Cohen’s class of distributions, the proposed method can extract high resolution estimates of both the envelope and the phase of the oscillations resulting in high resolution comodulograms. Unlike Hilbert transform based methods, the proposed method does not require bandpass filtering of the signals, and thus overcomes the problems of misidentification of the frequencies due to the selection of the bandwidth and the tran- sition band of the bandpass filter. While analytic wavelet transform based PAC methods address some of the shortcomings of Hilbert transform based PAC estimates, they are highly dependent on the selected wavelet parameters and the user provided range of frequencies. Unlike wavelet based PAC measures, the proposed approach does not depend on any design parameters and of- fers uniform frequency resolution thus the estimated PAC comodulograms do not change with the range of frequencies considered. Second, the proposed PAC measure is more robust to varying 44 signal conditions such as noise level, signal length, coupling strength and the separation between the coupled frequencies. In all of these cases, tf-MVL is more accurate at estimating the true cou- pling value. Finally, the proposed measure has been shown to be effective in detecting significantly higher alpha-gamma PAC for ERN compared to CRN. This finding is consistent with the results from previous studies [86, 91, 92], where an error related MEG study showed significant alpha gamma coupling. Dynamic coupling between alpha phase (8–13 Hz) and slower gamma ampli- tude (< 40Hz) has also been reported during visual processing in prior research [77, 92, 93]. It is hypothesized that during a visual task, the ongoing gamma-band activity in the visual cortex becomes temporally segmented by different phases of alpha-band activity [77, 92]. Moreover, the coupling strength was found to be higher for ERN compared to CRN, which is in line with the finding reported by Cohen et al. [94]. This can be explained by the fact that large-scale func- tional integration across different frequency bands supports flexible behavior adaption to improve the performance after an error [95]. The difference between ERN and CRN was most prominent at FCz, which is centered around the medial frontal cortex. This region was reported to be active during error processing and negative feedback [94] as it is hypothesized that error initiates the me- dial frontal based top-down control mechanisms to improve the performance [96, 97]. Botvinick et al. also reported that medial frontal cortex could detect error and conflicts and send signals to the lateral prefrontal cortex that behavioral adaptation is needed [98]. Thus, our findings are con- sistent with previous literature implicating higher alpha-gamma PAC in the medial frontal cortex and relating this with error-related performance. 45 2.6 Conclusion In this chapter, we evaluated the performance of a recently introduced time-frequency PAC method for estimating cross-frequency coupling. The proposed method offers higher resolution compared to existing Hilbert-based and wavelet-based methods and offers a framework to detect cross- frequency synchronization. Some of the limitations of the proposed time-frequency based am- plitude and phase estimation are as follows. First, the proposed method is based on a nonlinear quadratic transform of the signal as opposed to the Hilbert or wavelet transform based methods which are linear. As such, the computational complexity of the proposed method is O(N 2 log N) compared to Hilbert transform which has a complexity of O(N log N). Recently, it has been shown that quadratic time-frequency distributions can be approximated using linear transforms such as the short-time Fourier transform, thus obtaining computational complexity close to O(N log N) [99]. It is also worth mentioning that this computational complexity comes with the benefit that the proposed method does not depend on the choice of any parameters. The only parameter that can be tuned in the computation of RID-Rihaczek distribution is σ which controls the trade-off between resolution and cross-terms. In this chapter, σ = 0.001 was selected for all simulated and real data examples as it was observed that σ did not have a significant effect on the PAC values. In comparison, the Hilbert transform requires the choice of the filter parameters and the wavelet transform requires the choice of the parameters, β and γ, which greatly affect the amplitude and phase estimates. Second, the proposed method is more memory intensive compared to linear meth- ods as a complete time-frequency surface needs to be stored for extracting the amplitude and phase components at frequencies of interest. Due to this increased memory requirement and high com- putationally complexity, the proposed method can be applied to long epochs of neurophysiological recordings through a sliding window approach. Third, the proposed method similar to existing PAC 46 metrics is sensitive to data length. Simulation results indicate that the accuracy of the PAC value improves after 1000 samples. This is very similar to the behavior of the Hilbert transform based method. As such, the proposed method does not require more data samples than Hilbert transform to obtain accurate PAC estimates. Finally, the method presented in this chapter focuses on PAC within a signal, thus computing univariate PAC. However, most neurophysiological recordings in- volve data from multiple channels. Moreover, as indicated by de Cheveigne et al., a recording at a single electrode does not necessarily correspond to a single brain dynamic [100]. Therefore, a multivariate analysis is necessary to determine the cross-frequency interactions across channels as well as to separate the PAC within a channel from the nearby neural activity. Cohen [57] suggested a generalized eigenvalue decomposition based approach to address the latter. However, a complete bivariate PAC functional connectivity network analysis similar to bivariate PLV based connectivity networks [79] is still missing. The proposed univariate tf-MVL method can be easily extended to the bivariate case by computing PAC between all possible electrode and frequency pairs. This would result in a complete representation of the multi-frequency functional connectivity networks of the brain [101, 102]. The proposed tf-MVL measure can also be extended to consider different modes of coupling such as a biphasic coupling. 47 Chapter 3 Multivariate Analysis of Phase-Amplitude Coupling using Tensor Robust PCA Note: Most of the parts of this chapter have been published in the following article: ”Munia, Tamanna TK, and Selin Aviyente. ”Multivariate Analysis of Bivariate Phase-Amplitude Coupling in EEG Data Using Tensor Robust PCA.” IEEE Transactions on Neural Systems and Rehabilitation Engineering 29 (2021): 1268-1279.” 3.1 Introduction Neuronal oscillations and the interactions between them are believed to play a key role in under- standing the cognitive function of the brain [70, 72]. The interaction and coordination between oscillations at different frequencies is defined as cross-frequency coupling (CFC) [13, 14]. CFC plays a fundamental role in large scale neuronal encoding and communication by providing tempo- ral and spatial dynamics necessary for routing information through brain networks [13, 16]. CFC is an umbrella term to define a variety of phenomena like phase-phase coupling (PPC), amplitude- amplitude coupling (AAC), and phase-amplitude coupling (PAC). One of the most studied forms of CFC is PAC, which quantifies the interplay between the phase of a slower oscillation and the envelope of a faster oscillation [14]. PAC has been reported to 48 play a critical role in the execution of cognitive functions [103], attention selection [20], long-term memory processing [104], and sensory signal detection [52]. Although a variety of metrics have been developed to quantify PAC, most of these measures are bivariate in nature, i.e., they quantify local PAC observed between different frequencies of the same signal (within-channel PAC) or two different signals (cross-channel PAC). With the availability of multi-channel EEG data, there is a growing need for methods that can quantify neuronal couplings across the whole brain and across different frequency bands. Recently, several methods have been introduced to quantify multivariate phase amplitude cou- pling [75, 105, 106]. The existing methods are based on matrix factorization and as a result, they model multivariate PAC for a pre-defined pair of frequencies rather than capturing the variation of PAC across space and frequency, simultaneously. For example, phase coupling estimation (PCE) is limited to quantifying CFC between one high frequency and N low frequency signals. As such, PCE cannot capture the whole brain CFC dynamics. [105]. Generalized eigendecomposition (ged- CFC), on the other hand, is based on generalized eigendecomposition of multi-channel covariance matrices [75]. This is a hypothesis-driven framework and requires a priori knowledge about the low and high frequency bands that are coupled with each other. Thus, even though it is a multi- variate PAC analysis method, it is limited to two pre-determined frequency bands. In recent work [106], we have introduced a tensor based framework using PARAFAC decomposition to identify multivariate PAC. While this method can capture the multi-way nature of PAC across channels and frequency bands, it is highly dependent on model order selection. Moreover, it identifies the phase and amplitude providing channels separately rather than directly determining the coupled chan- nel pairs. Thus, this method requires additional significance testing of all possible combinations of phase and amplitude providing channels identified through the outer product of the extracted factors. 49 In this experiment, we first compute bivariate PAC values across all phase and amplitude pro- viding channels, frequency bands, and subjects and then represent these pairwise PAC values as a multi-way tensor. We propose a Higher order Robust Principal Component Analysis (HoRPCA) based multivariate framework to capture the spatial and spectral variation of the event-related PAC. HoRPCA decomposes an input tensor into low-rank and sparse parts and has been widely used in a variety of outlier detection applications such as distinguishing sparse event-related EEG from the background non-task related brain activity [107, 108, 109]. It has been previously reported that different behavioral tasks evoke distinct patterns of PAC across the cortex [12]. For example, enhanced PAC in the temporal cortex region is correlated with working memory maintenance tasks [110] whereas an increase in PAC in the medial frontal cortex is related to error processing tasks [111]. During an event-related study, as all subjects perform the same task and respond to the same stimulus, we assume a relatively similar PAC network structure exists for all subjects. A similar assumption was previously made about within frequency synchronization networks [109, 112]. Based on this assumption for event-related activity, HoRPCA based multivariate PAC approach captures the background PAC network in the low-rank part of the tensor while the sparse tensor captures the PAC activity associated with the event of interest. We evaluate the performance of HoRPCA based multivariate PAC approach on synthetic EEG data with varying signal parameters such as noise, subject variability, and volume conduction and compare it with existing multivariate PAC approaches. Finally, we apply the proposed method to EEG recordings collected during a cognitive control study to identify the spatial and spectral dynamics of pairwise PAC associated with different response types and to differentiate between the response types using multivariate analysis of PAC networks. 50 3.2 Background 3.2.1 Tensor Notation A multidimensional array with N modes given as X ∈ RI1 ×I2 ×...×IN is called a tensor, where xi1 ,i2 ,..iN denotes the (i1 , i2 , ..iN )th element. Column vectors obtained by fixing all indices of the tensor except the one that corresponds to nth mode are called mode-n fibers. The process of rearranging the tensor elements into a matrix is defined as unfolding or matricization. The mode-n unfolding of the tensor X , which is denoted by X(n) , can be obtained by arranging the mode-n fibers as the columns of the resulting matrix. The vectorization of X is defined as vec(X ). 3.3 Methods 3.3.1 Time-Frequency PAC (t-f PAC) PAC for a given signal is defined as the modulation between the amplitude A fa (t) of the high fre- quency ( fa ) with the phase φ f p (t) of the low frequency ( f p ) oscillations. The first step in computing PAC between two neuronal oscillations, x(t) and y(t) (where the amplitude of x(t) is coupled with the phase of y(t)) is to compute the RID-Rihaczek distributions Cx (t, f ) and Cy (t, f ) following (1.7). The amplitude component of x(t) at the desired frequency fa is extracted from the frequency f a2 f R constrained time marginal of Cx (t, f ) using (2.1) as Axa (t) = Cx (t, f )d f , where fa1 and fa2 refer f a1 to the bandwidth around the high frequency of interest, fa . The phase component of y(t) at the desired frequency f p is computed from the complex RID-Rihaczek distribution Cy (t, f ) using (2.2)   fp Cy (t, f p ) as φy (t) = arg |Cy (t, f p )| . After quantifying the amplitude and phase components, bivariate PAC between x(t) and y(t) can be computed using the amplitude normalized modulation index (MI) proposed by Özkurt and 51 Schnitzler [64] as follows: fp f jφy,k (t) ∑K a k=1 Ax,k (t)e 1 MIx,y ( f p , fa ,t) = √ q , (3.1) K K fa 2 ∑k=1 Ax,k (t) f a where K is the number of trials, Ax,k (t) is the extracted amplitude component from x(t) at high f frequency fa for the kth trial and φy,kp (t) is the extracted phase component at low frequency f p from y(t) for the kth trial. This metric quantifies PAC as a function of time, low frequency ( f p ) and high frequency ( fa ) and is between 0 and 1. As time-frequency based phase and amplitude component estimates were used to compute MIx,y ( f p , fa ,t), for the remainder of the chapter, this PAC index will be referred to as t-f PAC to differentiate it from conventional PAC metrics. MATLAB im- plementation of this metric is given in [113]. From this t-f PAC, we compute time averaged PAC ∑t∈T MIx,y ( f p , fa ,t) as: MIx,y ( f p , fa ) = |T | , where T is the time window of interest and |T | is the number of time points in that interval. This time averaged PAC, which is constructed by averaging the multi-dimensional t-f PAC estimate within a time window of interest, is used for constructing the PAC networks. As such, although t-f PAC itself is a time-varying metric, the PAC networks are not time-varying. 3.3.2 Statistical Significance Testing The significance of pairwise PAC values can be determined by surrogate data testing. The surrogate data is generated by following the block swapping procedure suggested in [114]. In this approach, a surrogate amplitude time series was generated by splitting the amplitude time series at a random point and swapping the two obtained time series. This surrogate amplitude time series was then associated with the original phase time series to compute PAC values. This procedure was repeated 52 100 times to generate 100 random time averaged surrogate PAC values for each frequency pair. Using these values, a threshold value, T hi, j ( f p , fa ) is obtained at the 95% confidence interval for each low frequency ( f p ) - high frequency ( fa ) combination between phase channel i and amplitude channel j. Time averaged PAC values that surpass this threshold value are considered significant for that low-high frequency pair and used for constructing the PAC network. We also obtain a threshold T hi, j (Fp , Fa ), for all 1 ≤ i, j ≤ N channel pairs and low and high frequency band pairs by averaging all the threshold values within those frequency bands, i.e., ∑ f p ∈Fp ∑ fa ∈Fa T hi, j ( f p , fa ) T hi, j (Fp , Fa ) = |Fp ||Fa | , where Fp is the low frequency band of interest, Fa is the high frequency band of interest, |Fp | is the number of frequencies in the low frequency band Fp and |Fa | is the number of frequencies in the high frequency band Fa . 3.3.3 PAC Tensor Construction Across Frequency Bands and Subjects For a channel and frequency band pair, we first compute the pairwise time averaged PAC value MIi, j ( f p , fa ) for all frequencies of interest, f p and fa and compare it with the threshold value T hi, j ( f p , fa ). MIi, j ( f p , fa ) values that surpass the corresponding threshold values are considered significant and used for computing the average PAC for that channel pair. After quantifying the pairwise time averaged PAC for all EEG channels and across all frequencies, a N × N weighted and F ,Fa directed connectivity network, AFp ,Fa , is constructed where, Ai jp = MIi, j (Fp , Fa ), for all 1 ≤ i, j ≤ ∑ f p ∈Fp ∑ fa ∈Fa MI ∗ i, j ( f p , fa ) N. Here, MIi, j (Fp , Fa ) = |Fp ||Fa | , where Fp is the low frequency band of interest, Fa is the high frequency band of interest, |Fp | is the number of frequencies in the low frequency band Fp , |Fa | is the number of frequencies in the high frequency band of interest Fa and    MI i, j ( f p , fa ), if MI i, j ( f p , fa ) > T hi, j ( f p , fa ),  ∗ MI i, j ( f p , fa ) =   0,  otherwise. 53 In this experiment, we focus on four different low frequency bands (Fp s), delta (1-3 Hz), theta (4-7 Hz), alpha (8-12 Hz), and beta (13-30 Hz) that may modulate the amplitude of the high frequency (Fa ) band gamma (31-100 Hz). After constructing these adjacency matrices for all fre- quency bands of interest, the goal of multivariate PAC analysis is to detect the spatial and spectral localization of the coupled channel pairs. Given AFp ,Fa , a 4-way tensor A ∈ RN×N×M×S is constructed where the first mode corresponds to N phase providing channels, second mode corresponds to N amplitude providing channels, third mode corresponds to the M different low-high frequency band pairs and the fourth mode of the tensor is included to capture the variation of PAC values across S subjects. 3.3.4 PARAFAC based Multivariate PAC We compare our proposed approach with PARAFAC based multivariate analysis. PARAFAC is a generalization of PCA for multi-way data and has been commonly used to extract information from multi-channel EEG data [115, 116, 117]. For PARAFAC based multivariate PAC, following [106], a 4-way PARAFAC decomposition is used to express A in terms of its factors across each mode. As A is non-negative, a non-negativity constraint was imposed on the PARAFAC factors as follows: R Ai jkl = ∑ λr air b jr ckr dlr , (3.2) r=1 where, air ≥ 0, b jr ≥ 0, ckr ≥ 0 and dlr ≥ 0 are the elements of the loading vectors across each mode. The first mode factors ar ∈ RN×1 provide the spatial loading for the different phase- providing channels whereas br ∈ RN×1 provide the spatial loading for the different amplitude- providing channels. The third mode factors, cr ∈ RK×1 correspond to the low frequency bands that modulate the high frequency band. The factors for the fourth mode dr ∈ RS×1 contain the signature 54 for each subject. As described in [106], the peaks in the factors of the first two modes allow us to determine the spatial locations of the phase- and amplitude-providing channels while the peaks in the factors of third mode provide the coupled low-high frequency band pairs. In this analysis, the rank, R was determined using DIFFerence in FIT (DIFFIT) method proposed in [118] as it was reported to be more suitable for offline analysis [119]. 3.3.5 HoRPCA based Multivariate PAC For matrix or two-way data, the limitations of PCA associated with outliers and non-Gaussian errors have been addressed using robust PCA (RPCA) [120]. In this method, a given matrix is decomposed into a low-rank and a sparse matrix. The extension of RPCA to higher-order data, i.e., tensors, has been proposed recently by Goldfarb and Qin [121]. In HoRPCA, the nuclear rank of a matrix is altered by the Tucker rank (Trank) of a tensor. To generalize RPCA to tensors, Trank and l0 norms are replaced with their convex surrogates CTrank and l1 norm. Given a tensor, X , the corresponding optimization problem is given as: minL ,S CTrank(L ) + λ kS k1 subject to L + S = X , (3.3) where L is the low-rank part of the tensor and S is the sparse part of the tensor and ||S ||1 := ||vec(S )||1 . In this experiment, we propose to employ HoRPCA [122, 123, 124] to capture the dynamics of multivariate PAC. HoRPCA is applied to A to obtain the estimates of low-rank component L and sparse component S . The low-rank tensor L accounts for the background or task-independent PAC network whereas the sparse component S represents the task-relevant PAC network. For computing HoRPCA, we considered the Singleton model proposed in [121], which evalu- 55 ates the nuclear norm of the tensor as the sum of the nuclear norms of the mode-n unfolding of the tensor, so CTrank(L ) := ∑4i=1 ||L(i) ||∗ . HoRPCA with the Singleton low-rank tensor model can be written as the following convex optimization problem: 4 minL ,S ∑ ||L(i)||∗+λ ||S ||1 subject to L + S = A . (3.4) i=1 The low-rank component L and sparse component S can be extracted from A by solving (3.4) using an alternating direction augmented Lagrangian (ADAL) method [121]. To apply HoR- PCA algorithm, we need to tune the regularization parameter λ , which controls the trade-off be- tween the sparse and low-rank parts of the tensor A . If λ is small, more data points will be allowed to move from L to S , and as λ increases fewer data points will be considered anomalous. For this study, following [121], λ was set as λ ≡ √1 where, Imax = max(N, N, M, S). Imax For multivariate PAC detection, we are interested in the sparse component S , as the nonzero elements of this tensor will provide us the phase and amplitude providing channels. Each row of S (:, :, m, l) will provide the phase providing channels while each column of S (:, :, m, l) will provide the amplitude providing channels for the mth frequency band and lth subject. Figure 3.1 il- lustrates the flowchart for the proposed HoRPCA based multivariate PAC detection method. A toy example detailing how the proposed HoRPCA based multivariate t-f PAC measure can detect cou- pled channel pairs and the frequency band pairs, simultaneously, is provided in the supplementary material. 56 Figure 3.1: Illustration of the proposed HoRPCA based multivariate t-f PAC measure. (a) Generation of the high-order tensor A by concatenating the N × N PAC networks across M various frequency bands and S subjects. (b) Low-rank (L ) and sparse tensor (S ) decomposition of A through HoRPCA. 3.4 Description of Datasets and Experiments The proposed multivariate PAC approach was validated first on multiple synthesized data sets, and then on EEG data collected during a cognitive control study. 3.4.1 Synthesized Multi-channel Data To generate synthesized EEG data, time series data were created at 2, 004 dipole locations in the brain. These dipole locations were based on the standard MRI brain. EEG data were created by considering the fact that the spectrum of EEG data follow the power law, i.e., higher the fre- quency, weaker the amplitude (P( f ) ∼ ( f1β )). The rate of decay of the amplitude is defined by the parameter β with β = 0, 1, and 2 indicating white noise, pink noise, and Brownian (red) noise, respectively [125]. The frequency spectrum of the human brain has been reported to be similar to Brownian noise [81, 82], hence EEG data was generated from Brownian noise. Time series data from Brownian noise was generated at each of the 2, 004 uncorrelated dipole locations in the brain 57 at a sampling frequency of Fs Hz. Cross-voxel correlations were imposed across all the dipoles by generating a random dipole-to-dipole correlation matrix. The new data matrix was constructed as Y = X T V D1/2 , where V is the eigenvector, D is the eigenvalue, and X is the data matrix. To project each dipole location to the scalp EEG locations, the forward model was generated by following the openmeeg [126] algorithm implemented in Brainstorm [127]. To introduce phase-amplitude coupling in the simulated data, first the dipole locations corre- sponding to the high frequency amplitude component and low frequency phase component were chosen. Brownian noise signal xB (t) was generated with a sampling frequency of Fs Hz. Low frequency phase signal x f p (t) was generated by bandpass filtering xB (t) at the phase providing frequency f p with a bandwidth of 2 Hz ( f p ± 1Hz). To create the amplitude signal, xB (t) was bandpass filtered at the amplitude providing frequency fa with a bandwidth of 10 Hz ( fa ± 5Hz) to ensure that the high frequency activity is broadband. PAC was generated using the procedure de- scribed by Kramer and Eden [84]. The time locations of relative maxima of the phase signal x f p (t) were detected. At each maxima, a DC shifted Hanning window with a duration of 42 ms, was multiplied with the amplitude time series. The amplitude signal x fa (t) with monophasic coupling was generated at high frequency fa by multiplying the Hanning window with the filtered ampli- tude time series centered at the relative maxima (peaks) of the phase signal x f p (t). The coupling intensity is controlled by multiplying the Hanning window itself with a multiplier I, where I = 1 indicates full coupling and I = 0 indicates no coupling. Finally, the time series of the chosen dipole locations were replaced with x f p (t) for the low frequency and x fa (t) for the high frequency signal. The whole procedure was repeated d times to generate d distinctly coupled neuronal oscillations between 2d dipole locations (d phase providing dipoles and d amplitude providing dipoles). 58 3.4.2 EEG Data The EEG dataset used to assess the proposed PAC measure was collected during a cognitive control-related error processing study. The experiment conducted was a letter version of the speeded-reaction Flanker task [128]. The study was designed following the experimental protocol and guidelines approved by the Institutional Review Board (IRB) of the Michigan State University. The data acquisition was performed following the guidelines and regulations established by this protocol. Written and informed consent was collected from each participant before data collec- tion. A string of five letters, which could be congruent (e.g., SSSSS) or incongruent (e.g., SSTSS), was exhibited in front of each participant for each trial. The participants have to select the center letter with a standard mouse with respect to the Flanker letters. The trials began with a 35ms of flanking stimuli (e.g., SS SS). Then, the target stimuli were embedded in the center of the flankers (e.g., SSSSS/SSTSS) and remained for 100 ms (total presentation time is 135ms) followed by an inter-trial interval of varying duration ranging from 1200 to 1700 ms. Continuous EEG responses were recorded using the ActiveTwo system (BioSemi, Amsterdam, The Netherlands) by placing the 64 Ag–AgCl channels following the international 10/20 system. All EEG signals were sampled at 512 Hz. The trials with artifacts were removed, and volume conduction was minimized using the Current Source Density (CSD) Toolbox [90]. The artifact removed data were considered for analysis in this experiment. The trials corresponding to the error and correct responses were separated from each other for analysis. Each trial was one second long. A total of 20 participants were considered for the ongoing study. The inclusion criteria were that the number of trials for error response should be at least 20 or higher. Since the total number of correct/error response trials varies for different participants, the number of trials considered for both responses were kept the same for a fair comparison. The correct response trials were randomly 59 sub-sampled from the whole set of correct responses and the whole procedure was repeated 10 times to select the correct response trials. 3.4.3 Synthesized Data Experiments The proposed HoRPCA based multivariate PAC method is first tested on the synthesized data set described in Section 3.4.1. Experiments were conducted to evaluate the robustness of the proposed method against spurious coupling due to noise, subject variability and linear mixing with respect to existing approaches. In these experiments, HoRPCA based method was compared to three other methods, i.e., simple averaging, matrix factorization based gedCFC [75] and PARAFAC based method [106]. For simple averaging method, we first average the delta-gamma, theta-gamma, alpha-gamma and beta-gamma PAC networks for each subject. We also average the thresholds for delta-gamma, theta-gamma, alpha-gamma and beta-gamma networks obtained from the statistical significance testing described in Section 3.3.2 for each subject. To detect the coupled channels, we compare the averaged PAC network with the thresholds for each electrode pair. PAC values for gedCFC were computed following the codes provided in [75]. As gedCFC computes PAC between two pre-defined low and high-frequency bands, we compute gedCFC based PAC individually for each frequency band pair combination. PARAFAC based analysis is conducted as described in section 3.3.4. For both gedCFC and PARAFAC based methods, to differentiate the coupled chan- nel pairs from background PAC values, the outputs of these two methods are compared to the average threshold PAC value (T hi, j (Fp , Fa )) for that channel pair obtained through the statistical significance testing described in Section 3.3.2. In this section, we conducted two experiments. Experiment 1 focuses on whether the methods can differentiate true coupling from no coupling whereas experiment 2 focuses on the accuracy of the detected channel pairs. The performance evaluation metrics and the description of the experi- 60 ments are as follows. Performance evaluation metrics: The performance metrics are: TP 1. Precision = T P+FP , TP 2. Recall/sensitivity = T P+FN , 2×precision×recall 3. F-measure = precision+recall , √ 4. G-mean = sensitivity × speci f icity, TN where, Specificity = T N+FP , where, • TP (true positive): coupled pairs detected as coupled. • FP (false positive): uncoupled pairs detected as coupled. • FN (false negative): coupled pairs detected as uncoupled. • TN (true negative): uncoupled pairs detected as uncoupled. Experiment 1: Evaluation of detection power of different approaches In this experiment, we are interested in quantifying the ability of the different methods to differen- tiate between multi-channel data with and without PAC for different noise levels. This simulation is conducted to ensure that our method is robust against external noise, e.g., measurement or bi- ological noise, and is able to differentiate true coupling from uncoupled oscillatory activity. Two multi-channel EEG datasets were generated. First dataset was generated with theta-gamma cou- pling (I = 0.7) between any two randomly selected dipole locations. Signals for amplitude and phase providing dipoles were generated following the procedure described in Section IV-A. For all 61 other dipoles, the signals were time series data generated from Brownian noise. This procedure was repeated 20 times to generate a dataset with 20 subjects with different dipole locations and this dataset was referred to as ’dataset with coupling’. For the second dataset, low and high frequency signals were generated at the same dipole locations as the first dataset but the coupling intensity was set to zero (I = 0) to generate no-coupling condition. This procedure was also repeated 20 times to generate a dataset with 20 subjects and this dataset was referred to as ’dataset without coupling’. The performance of each approach was determined by evaluating whether each method detected any coupled pairs for each subject or not. Based on this evaluation, detection accuracy in terms of precision, recall, F-measure and G-mean is quantified across subjects for varying signal to noise ratio (SNR) levels (-6dB to 6dB). We repeat the whole experiment 20 times to quantify the variability of each method. Experiment 2: Performance evaluation with varying signal parameters We next conduct experiments assessing the accuracy of the detected coupled channel pairs with respect to various signal parameters, including noise, variability in dipole locations and number of coupled channel pairs. As gedCFC can not detect multiple channel pairs across multiple frequen- cies, it was not included in these comparisons. For evaluating the performance of the remaining three methods (HoRPCA, PARAFAC and averaging), for a fair comparison, the number of channel pairs detected by each method was fixed as the actual number of coupled pairs in the data. Robustness to noise: To evaluate the robustness of the proposed method, we generated synthe- sized data with additive white Gaussian noise of varying SNR from -6dB to 6dB added to all dipole locations. Assessing the robustness to noise is important as real EEG signals are not simple linear combinations of dipole sources and noise can be generated by both neural and external measure- ment systems. Eight coupled pairs (5 theta-gamma and 3 alpha-gamma) were generated for each subject between 16 dipole locations (8 phase providing and 8 amplitude providing dipoles) fol- 62 lowing the steps described in Section 3.4.1. The projections of these dipole locations are shown in Figure 3.2. This procedure was repeated ten times to generate a dataset with ten subjects for each SNR level. The whole experiment was repeated 20 times to quantify the variability of each method. Figure 3.2: The projections of the dipole locations to the scalp for the phase and amplitude components of the synthesized data for the performance evaluation experiments with varying signal parameters. (a) Theta band phase providing dipoles; (b) Gamma band amplitude providing dipoles coupled with the theta phase providing dipoles; (c) Alpha band phase providing dipoles; (d) Gamma band amplitude providing dipoles coupled with the alpha phase providing dipoles. Robustness to subject variability: For multi-subject analysis, methods such as our method which tries to identify coupled channels across subjects, variability across subjects is natural. In this experiment, we evaluate the performance of the proposed method by introducing variability in the location of the coupled channels to determine the robustness of the method against variability across subjects. First, eight multivariate coupled pairs of oscillations (5 theta-gamma and 3 alpha- gamma) were generated for each subject between 16 locations (8 phase providing and 8 amplitude providing) following the steps described in Section 3.4.1 and for the dipoles illustrated in Figure 3.2. Using this procedure, a dataset was generated with ten subjects. SNR was fixed at 6dB for all subjects. Variability was introduced by changing 5% to 25% of the dipole locations with a step size of 5% for different datasets. We repeated the whole experiment 20 times to quantify the variability of each method. Robustness to volume conduction: The performance of the proposed measure was also evalu- 63 ated with increasing number of coupled channel pairs to determine the robustness of the proposed method against volume conduction. Volume conduction is very commonly encountered in EEG and refers to the phenomenon where nearby electrodes are highly correlated with each other due to the conductivity of the scalp [129]. This is important as the number of coupled dipoles increases there will be more interference across channels which may make it harder to detect the truly cou- pled channels. Synthesized datasets (10 subjects in each dataset with the same number of coupled pairs) were generated by varying the number of coupled channel pairs from 4 to 20 following the steps described in Section 3.4.1. SNR was fixed at 6dB for all subjects. 75% of the coupled channel pairs were theta-gamma coupled while 25% were alpha-gamma coupled. The whole experiment was repeated 20 times to quantify the variability of each method. 3.4.4 EEG Data Experiments Detection of PAC networks Using the average PAC within a frequency band and time window of interest (MI) defined in section III-C, pairwise PAC values were computed across all 64 electrodes for delta-gamma, theta- gamma, alpha-gamma and beta-gamma frequency band combinations for each subject and re- sponse type, error and correct. As previous studies indicate increased synchronization associated with the ERN in the time window 0-100 ms [79], the connectivity matrices were constructed by averaging the time-varying PAC within this time window. This resulted in Aerror ∈ R64×64×4×20 and Acorrect ∈ R64×64×4×20 . HoRPCA is applied on both Aerror ∈ R64×64×4×20 and Acorrect ∈ R64×64×4×20 to obtain the estimates of low-rank components Lerror , Lcorrect and sparse components Serror , Scorrect . The sparse components indicate that the theta-gamma frequency band pairs have the highest PAC for both EEG error and correct responses. The set of coupled channels for error and correct responses 64 for theta band across all subjects was determined for all 1 ≤ i, j ≤ N as: 20 Error response: Xer (i, j) = ∑ 1R+ [Serror (i, j, 2, l)], (3.5) l=1 20 Correct response: Xcr (i, j) = ∑ 1R+ [Scorrect (i, j, 2, l)], (3.6) l=1   1, if Serror (i, j, 2, l) > 0,   where, 1R+ [Serror (i, j, 2, l)] = and   0, otherwise,    1, if Scorrect (i, j, 2, l) > 0.   1R+ [Scorrect (i, j, 2, l)] =   0, otherwise.  The detected channel combinations for error (Xer ) and correct (Xcr ) responses are referred to as PAC networks for error and correct responses, respectively. Classification based on detected PAC networks A classification experiment was conducted to determine if the detected multivariate PAC networks can be used as features to discriminate between error and correct responses. As the proposed method uses the sparse part of the tensor for classification and since this sparse part is data depen- dent, the classification was performed by following the sparse representation classification (SRC) approach described in [130, 131]. We first generate training and test samples from Aerror ∈ R64×64×4×20 and Acorrect ∈ R64×64×4×20 using 5-fold cross-validation. For each validation, we take I subjects (both error and correct re- sponses) to generate a test set Atest ∈ R64×64×4×2I = [Aerrtest , Acrrtest ]. We then compute HoRPCA for the remaining J = 20 − I subjects in the training set separately for error Aerrtrain ∈ R64×64×4×J and correct responses Acrrtrain ∈ R64×64×4×J to obtain two separate training feature matrices cor- 65 responding to error and correct responses. The feature matrices are the fourth mode slices of the sparse tensors Serrtrain ∈ R64×64×4×J and Scrrtrain ∈ R64×64×4×J corresponding to each sub- ject j ∈ {1, 2, ..., J} in the training set. Next, for each subject and response type in the test set (Yk = Atest (:, :, :, k)), we concatenate the error and correct tensors with the PAC tensor corre- sponding to each test sample (Yk ) separately, i.e., we end up with two tensors [Aerrtrain , Yk ] ∈ R64×64×4×(J+1) and [Acrrtrain , Yk ] ∈ R64×64×4×(J+1) . We perform HoRPCA on each of these ten- sors and extract the sparse matrix corresponding to the test sample. This results in two feature matrices for our test sample, one based on projection with respect to the error training samples and the other based on projection with respect to correct training samples. For each of these feature matrices, we compute the distances between the feature matrix from the test sample and the feature matrices obtained from the training data. The test sample is assigned to the class to which it has the smallest distance, i.e., Nearest Neighbor (NN) classification. We repeat this procedure across all test samples and multiple 5-folds and report the average accuracy. The overall procedure is summarized in Algorithm 1. We compare the classification performance of HoRPCA with simple averaging, gedCFC and PARAFAC based methods. As simple averaging and gedCFC methods are not data dependent, to classify the error and correct responses for these two methods, we employ simple NN classifier without the SRC based algorithm. For simple averaging method, the classification was performed using 5-fold cross-validation with an NN classifier. The feature matrices in this case are the 64×64 PAC networks averaged across frequency bands corresponding to each subject and response type. The classification for gedCFC was also performed by using 5-fold cross-validation with NN clas- sifier. In this case, the features are the multivariate PAC values for the theta-gamma band extracted by gedCFC. For PARAFAC based method, as it is also data dependent similar to HoRPCA, the classification was performed following the SRC based NN-classifier approach described in Algo- 66 Algorithm 1: Sparse Representation-based Classification (SRC) of PAC Networks Input : Training tensors: Aerrtrain ∈ R64×64×4×J and Acrrtrain ∈ R64×64×4×J ; Test tensor: Atest ∈ R64×64×4×2I = [Aerrtest , Acrrtest ], where, Aerrtest ∈ R64×64×4×I , Acrrtest ∈ R64×64×4×I . Output: Label of test samples in Atest . 1 Perform HoRPCA on Aerrtrain and Acrrtrain using (3.4) to obtain training features Serrtrain and Scrrtrain . 2 for k = 1 to 2I do 3 Yk = Atest (:, :, :, k). 4 Using (3.4), perform HoRPCA on A ∗ errtrain = [Aerrtrain , Yk ], A ∗ crrtrain = [Acrrtrain , Yk ] to obtain S ∗ errtrain and S ∗ crrtrain . 5 Serr k test ← S ∗ errtrain (:, :, 2, J + 1) 6 Scrr k test ← S ∗ crrtrain (:, :, 2, J + 1). 7 for l = 1 to J do 8 l dk,err = ||Serrtrain (:, :, 2, l) − Serrk || test 2 9 dk,crr = ||Scrrtrain (:, :, 2, l) − Scrrtest ||2 l k 10 end 11 Label(Yk ) = argminerr,crr (dk,err l l , dk,err ) for all l. 12 end rithm 1. For all cases, the performance was evaluated in terms of precision, recall, F-measure and G-mean. 3.5 Results 3.5.1 Synthesized Data Results Results for Experiment 1 The evaluation performance of the detection power for the four methods for various SNR values is reported in Table 3.1. From Table 3.1, it can be seen that the proposed HoRPCA based PAC ap- proach outperforms gedCFC in all cases and can detect the presence of PAC even for highly noisy data (-6dB) with a F-measure close to 0.85. The detection power of PARAFAC based approach was also high and comparable to the proposed HoRPCA based approach. The high performance 67 exhibited by the HoRPCA and PARAFAC based methods indicates that the tensor based multi-way analysis takes full advantage of the multi-linear structure of the data, which improves the classifi- cation performance compared to matrix factorization-based methods like gedCFC. Averaging the networks has the worst performance (i.e., performance around the chance level) for all SNR levels indicating loss of information due to averaging. Table 3.1: Comparison of PAC detection power obtained with HoRPCA, gedCFC, PARAFAC and Averaging based methods for different SNR values. The table depicts the Mean±standard deviation of evaluation matrices for 20 repetition of each experiment. Here, Prec.: Precision, Rec.:Recall, F-meas.: F-measure. SNR Multivariate PAC with Averaging gedCFC based multivariate PAC PARAFAC based multivariate PAC HoRPCA based multivariate PAC F- G- F- G- F- G- F- G- (dB) Prec. Rec. Prec. Rec. Prec. Rec. Prec. Rec. meas. mean meas. mean meas. mean meas. mean 0.37± 0.38± 0.38± 0.36± 0.75± 0.78± 0.77± 0.76± 0.79± 0.83± 0.82± 0.81± 0.81± 0.89± 0.85± 0.84± -6dB 0.112 0.114 0.113 0.115 0.083 0.092 0.090 0.092 0.077 0.079 0.078 0.079 0.067 0.065 0.063 0.062 0.43± 0.48± 0.46± 0.43± 0.79± 0.83± 0.81± 0.81± 0.84± 0.88± 0.87± 0.86± 0.86± 0.94± 0.89± 0.89± -3dB 0.113 0.109 0.111 0.113 0.088 0.081 0.081 0.080 0.073 0.073 0.073 0.070 0.060 0.062 0.059 0.059 0.50± 0.50± 0.49± 0.50± 0.84± 0.88± 0.86± 0.86± 0.89± 0.94± 0.91± 0.91± 0.92± 0.97± 0.95± 0.94± 0dB 0.100 0.103 0.101 0.103 0.078 0.081 0.079 0.081 0.063 0.063 0.064 0.063 0.062 0.051 0.052 0.051 0.59± 0.56± 0.57± 0.58± 0.91± 0.95± 0.94± 0.93± 0.98± 0.95± 0.97± 0.97± 0.96± 0.97± 0.98± 0.98± 3dB 0.090 0.093 0.092 0.093 0.068 0.053 0.059 0.055 0.050 0.043 0.043 0.043 0.042 0.042 0.040 0.034 0.62± 0.56± 0.59± 0.60± 0.94± 0.96± 0.95± 0.94± 1.00± 1.00± 1.00± 1.00± 1.00± 1.00± 1.00± 1.00± 6dB 0.083 0.093 0.091 0.083 0.065 0.060 0.059 0.053 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 Results for Experiment 2 The performance of the three approaches for various SNR levels, variability in dipole locations and number of coupled channel pairs are depicted in Tables 3.2, 3.3 and 3.4, respectively. From Tables 3.2, 3.3 and 3.4, we can see that the proposed HoRPCA based method outperforms PARAFAC and the averaging based methods in all cases. HoRPCA detected the coupled channels correctly with more than 80% precision and recall even for high noise levels, e.g. -6dB, (Table 3.2). HoRPCA is also shown to be more robust against subject variability (Table 3.3). It is also less sensitive to the number of coupled channel pairs (Table 3.4). While PARAFAC based multivariate PAC detection performs well for low noise levels and less variability in the data, its performance deteriorates quickly as the uncertainty in the data increases. For example, PARAFAC is unable to detect the correct coupled channel pairs in the presence of variability across subjects. This can be explained by the fact that PARAFAC based method detects 68 the channel pairs which are common across all subjects. In all cases, averaging the networks yielded the worst performance and was found to be more sensitive to different signal parameters. This is due to the fact that averaging cannot capture variability across subjects and frequency bands. The standard deviations reported in these Tables reflect that the variability of the HoRPCA method is smaller than others across multiple realizations of the experiments. This indicates the stability of the proposed method. Table 3.2: Comparison of detected channel pair accuracy for various SNRs obtained with HoR- PCA, PARAFAC and Averaging based methods. The table depicts the Mean±standard deviation of evaluation matrices for 20 repetition of each experiment. Here, Prec.: Precision, Rec.:Recall, F-meas.: F-measure. SNR Multivariate PAC with Averaging PARAFAC based multivariate PAC HoRPCA based multivariate PAC F- G- F- G- F- G- (dB) Prec. Rec. Prec. Rec. Prec. Rec. meas. mean meas. mean meas. mean 0.38± 0.38± 0.38± 0.61± 0.70± 0.70± 0.70± 0.84± 0.82± 0.82± 0.82± 0.90± -6dB 0.103 0.103 0.103 0.081 0.093 0.093 0.093 0.072 0.071 0.071 0.071 0.054 0.48± 0.48± 0.0.48± 0.69± 0.79± 0.79± 0.79± 0.88± 0.84± 0.84± 0.84± 0.92± -3dB 0.099 0.099 0.099 0.076 0.056 0.056 0.056 0.047 0.046 0.046 0.046 0.024 0.55± 0.55± 0.55± 0.74± 0.84± 0.84± 0.84± 0.90± 0.90± 0.90± 0.90± 0.94± 0dB 0.093 0.093 0.093 0.057 0.056 0.056 0.056 0.039 0.031 0.031 0.031 0.019 0.61± 0.61± 0.61± 0.78± 0.92± 0.92± 0.92± 0.95± 0.95± 0.95± 0.95± 0.97± 3dB 0.076 0.076 0.076 0.047 0.054 0.054 0.054 0.028 0.028 0.028 0.028 0.015 0.66± 0.66± 0.66± 0.81± 0.96± 0.96± 0.96± 0.98± 0.99± 0.99± 0.99± 0.99± 6dB 0.086 0.086 0.086 0.054 0.018 0.018 0.018 0.009 0.004 0.004 0.004 0.002 3.5.2 EEG Data Results PAC network detection results The low-rank and sparse networks extracted for the theta band of error and correct responses averaged across subjects is shown in Figure 3.3. This figure shows that the low-rank parts of the tensor for both response types are similar, while the sparse parts illustrate the differences in the coupled brain regions for the two response types. This observation is also verified by statistical testing. A Wilcoxon Signed Rank Sum Test (p < 0.05) with Bonferroni correction is conducted to determine the statistical significance of the 69 Table 3.3: Comparison of detected channel pair accuracy for variability in channel locations obtained with HoRPCA, PARAFAC and Averaging based methods. The table depicts the Mean±standard deviation of evaluation matrices for 20 repetition of each experiment. Here, Prec.: Precision, Rec.:Recall, F-meas.: F-measure. Variability Multivariate PAC with Averaging PARAFAC based multivariate PAC HoRPCA based multivariate PAC F- G- F- G- F- G- (%) Prec. Rec. Prec. Rec. Prec. Rec. meas. mean meas. mean meas. mean 0.66± 0.66± 0.66± 0.81± 0.96± 0.96± 0.96± 0.98± 0.99± 0.99± 0.99± 0.99± 0% 0.086 0.086 0.086 0.054 0.018 0.018 0.018 0.009 0.004 0.004 0.004 0.002 0.61± 0.61± 0.61± 0.78± 0.90± 0.90± 0.90± 0.95± 0.97± 0.97± 0.97± 0.98± 5% 0.076 0.076 0.076 0.062 0.029 0.029 0.029 0.016 0.011 0.011 0.011 0.006 0.56± 0.56± 0.56± 0.75± 0.82± 0.82± 0.82± 0.90± 0.95± 0.95± 0.95± 0.97± 10% 0.084 0.084 0.084 0.063 0.038 0.038 0.038 0.023 0.026 0.026 0.026 0.013 0.44± 0.44± 0.44± 0.66± 0.77± 0.77± 0.77± 0.87± 0.92± 0.92± 0.92± 0.96± 15% 0.090 0.090 0.090 0.063 0.047 0.047 0.047 0.025 0.028 0.028 0.028 0.016 0.40± 0.40± 0.40± 0.63± 0.72± 0.72± 0.72± 0.84± 0.90± 0.90± 0.90± 0.95± 20% 0.097 0.097 0.097 0.065 0.055 0.055 0.055 0.034 0.032 0.032 0.032 0.017 0.34± 0.34± 0.34± 0.58± 0.68± 0.68± 0.68± 0.82± 0.86± 0.86± 0.86± 0.92± 25% 0.105 0.105 0.105 0.067 0.070 0.070 0.070 0.043 0.044 0.044 0.044 0.024 difference between ERN and CRN PAC networks across subjects. The statistical testing revealed no significant difference between the low rank networks but significant difference was observed between sparse networks. We constructed a p-value PAC network that shows the channel pairs with significant difference between the sparse networks corresponding to error and correct responses. The p-value PAC network is depicted in Figure 3.4. 70 Table 3.4: Comparison of detected channel pair accuracy for various number of coupled channel pairs obtained with HoRPCA, PARAFAC and Averaging based methods. The table depicts the Mean±standard deviation of evaluation matrices for 20 repetition of each experiment. Here, Prec.: Precision, Rec.:Recall, F-meas.: F-measure. Coupled Multivariate PAC with Averaging PARAFAC based multivariate PAC HoRPCA based multivariate PAC F- G- F- G- F- G- Pairs Prec. Rec. Prec. Rec. Prec. Rec. meas. mean meas. mean meas. mean 0.72± 0.72± 0.72± 0.85± 1.00± 1.00± 1.00± 1.00± 1.00± 1.00± 1.00± 1.00± 4 0.092 0.092 0.092 0.0062 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.66± 0.66± 0.66± 0.81± 0.96± 0.96± 0.96± 0.98± 0.99± 0.99± 0.99± 0.99± 8 0.086 0.086 0.086 0.054 0.018 0.018 0.018 0.009 0.004 0.004 0.004 0.002 0.56± 0.56± 0.56± 0.75± 0.90± 0.90± 0.90± 0.95± 0.96± 0.96± 0.96± 0.98± 12 0.087 0.087 0.087 0.059 0.041 0.041 0.041 0.022 0.020 0.020 0.020 0.010 0.49± 0.49± 0.49± 0.70± 0.85± 0.85± 0.85± 0.92± 0.92± 0.92± 0.92± 0.96± 16 0.092 0.092 0.092 0.065 0.053 0.053 0.053 0.039 0.027 0.027 0.027 0.014 0.47± 0.47± 0.47± 0.68± 0.81± 0.81± 0.81± 0.89± 0.88± 0.88± 0.88± 0.93± 20 0.103 0.103 0.103 0.070 0.076 0.076 0.076 0.053 0.056 0.056 0.056 0.030 Figure 3.3: Low-rank and sparse networks extracted for the error and correct responses averaged across subjects. (a) Average low-rank PAC network for error response; (b) Average low-rank PAC network for correct response; (c) Average sparse PAC network for error response; (d) Average sparse PAC network for correct response. 71 Figure 3.4: PAC network of p-values showing significant difference between error and correct responses. The color scale from white to black indicates the level of significance (Wilcoxon Signed Rank Sum Test with p < 0.05). The arrows originate from phase providing channel. From the extracted sparse networks, we detected the channel combinations for error (Xer ) and correct (Xcr ) responses as shown in Figure 3.5 (a) and (b), respectively. The colormap shows the number of subjects for which a channel combination was detected. The yellow pixels indicate channel pairs with coupling across the majority of subjects, while the light blue pixels correspond to channel pairs with coupling for only a few number of subjects. 72 Figure 3.5: The detected channel combinations across all subjects for (a) error response (Xer ); and (b) correct response (Xcr ) for theta frequency band. The channel pairs (phase-amplitude coupled channels) detected for a majority of subjects for error response are: FP1-POz, FP1-Pz, FP1-C4, F3-Oz, FT7-POz, FC1-P7, FC1-POz, FC1-Pz, 73 FC1-C4, FC1-P6, Fz-C4, Fz-P6, FCz-P7, FCz-POz, FCz-Pz, FCz-P6, C2-Pz, and C2-P6. Sim- ilarly, for correct response the channel pairs (phase-amplitude coupled channels) detected for a majority of subjects are: FP1-P7, FP1-P3, FP1-O2, F1-Pz, F5-PO4, FPz-P7, FPz-P3, FPz-P8, FPz-C6, FPz-TP8, AFz-PO4 and C6-Pz. The resulting theta-gamma band PAC networks for error and correct responses are shown in Figure 3.6(a) and (b) respectively. Figure 3.6: Theta-gamma PAC networks detected for (a) Error response and (b) Correct response. The color scale from white to black indicates the number of subjects for which a channel pair was detected. The arrows originate from phase providing channel. Classification experiment result Table 3.5 shows the comparison of classification performance for 50 repetitions of the 5-fold cross- validation in terms of precision, recall, F-measure and G-mean. From this Table, it can be seen that the proposed HoRPCA based multivariate method results in very high accuracy (F-measure and G-mean > 0.98) compared to averaging (F-measure and G-mean close to 0.70) and matrix factorization-based gedCFC (F-measure and G-mean close to 0.85). The discrimination power of PARAFAC was slightly lower than the proposed method. The increased discrimination power of the tensor based methods indicates that the tensor based multivariate t-f PAC approaches are more 74 effective at identifying the brain regions central to PAC. Table 3.5: Classification of EEG error and correct responses based multivariate PAC networks computed using averaging, gedCFC, PARAFAC and HoRPCA based methods. Method Precision Recall F-measure G-mean Averaging 0.685±0.243 0.726±0.230 0.678±0.199 0.701±0.163 gedCFC 0.801±0.178 0.915±0.137 0.833±0.115 0.842±0.110 PARAFAC 0.906±0.141 1.000±0.00 0.945±0.0879 0.946±0.0797 HoRPCA 0.975±0.061 1.000±0.00 0.986±0.033 0.984±0.038 3.5.3 Interrelationship between Multivariate PAC and hormone level We also investigated the efficacy of inter-areal PAC for characterizing cross-frequency coupling in event-related potential data related to cognitive control. We investigated the interrelationship between inter-areal PAC with the impact of anxiety in women by examining the interplay of cog- nitive control dysfunction and ovarian hormones with the quantified PAC value. In particular, we examined whether and how estradiol modulated error-related inter-areal theta-gamma PAC. As FCz-POz provided the highest inter-channel t-f PAC for error-related activity, we conducted this experiment for this channel pairs. Inter-areal PAC value was computed for 43 subjects across four visits. Data collected during the four visits were assigned to different groups based on the partici- pant’s estradiol level during the data collection. Four groups were referred to as Level 1 to Level 4, with Level 1 indicating the lowest and Level 4 indicating the highest estradiol value during the data collection. For each visit, HoRPCA was applied to the resulting tensors and channel pairs with high coupling values were subsequently detected to compute the associations between the error-related PAC measure and endogenous estradiol levels. To determine if the computed theta- gamma interchannel t-f PAC values are significantly different across estradiol levels, a Wilcoxon Signed Rank Sum test with α = 0.05 was conducted. Finally, we conducted a correlation analy- sis to examine the interrelationship between the measured theta-gamma interchannel t-f PAC with 75 estradiol values. Theta-gamma interchannel t-f PAC between FCz-POz channel pair across 15 trials for 43 par- ticipants are depicted in Figure 3.7 through box plots. This figure shows an increase in theta- gamma interchannel t-f PAC with the increase of hormone levels across visits. Figure 3.7: Box plot of the theta-gamma interchannel t-f PAC between FCz-POz channel pair across different hormone levels (N = 43). The central red lines indicate the median, and the bottom and top edges of the boxes indicate the 25th and 75th percentiles, respectively. The whiskers extend till the most extreme data points, which are not considered outliers, and the outliers are plotted individually by the ’+’ symbol. The Wilcoxon Signed Rank Sum test with α = 0.05 revealed some significant difference across estradiol levels in terms of theta-gamma coupling. The p-values corrected through Bonferroni correction are given in Table 3.6 for all hormone levels. As demonstrated in Table 3.6, the average theta-gamma interchannel t-f PAC increased significantly with an increase of estradiol level (p < 0.01 between Levels 1 and 4 and p < 0.05 between Levels 2 and 4). 76 Table 3.6: Comparison of theta-gamma interchannel t-f PAC across different estradiol levels. Bold values indicate significant values after Bonferroni correction. Here, * indicates significantly dif- ferent for α = 0.05 and ** indicates significantly different for α = 0.01. p values Level 1 Level 2 Level 3 Level 4 Level 1 1.000000 0.636116 0.068551 0.001076** Level 2 0.636116 1.000000 0.453425 0.026779* Level 3 0.068551 0.450934 1.000000 0.453425 Level 4 0.001076** 0.026779* 0.453425 1.000000 A scatter plot showing the distribution of the measured theta-gamma interchannel t-f PAC with estradiol values is given in Figure 3.8. The correlation between estradiol values and computed t-f PAC values was calculated by Pearson Correlation Matrix and is shown in Figure 3.8. As shown in Figure 3.9, correlation analysis showed that the theta-gamma t-f PAC values between FCz-POz channels were highly correlated (r = 0.43, p < 0.001) with increased estradiol levels. Figure 3.8: Scatter-plot depicting the relationship between the computed theta-gamma t-f PAC and estradiol values. 77 Figure 3.9: Pearson correlation analysis for theta-gamma t-f PAC and estradiol values. 3.6 Discussion In this chapter, a HoRPCA based approach was proposed to identify the spatial and spectral com- ponents of multivariate PAC across all channels, frequency bands, and subjects. As illustrated in Figure 3.5, HoRPCA can determine the coupled channel pairs directly from the row and column indices of the non-zero elements of the sparse tensor component for a given frequency band. This approach provides certain advantages over PARAFAC decomposition of the same tensor. First, with the PARAFAC model, order selection is an open problem. While the first factor across each mode can be used to identify the coupled channel pairs and frequency bands, these factors are only effective at capturing the largest variance in the data and are not necessarily robust to outliers and variations in the data such as the ones considered in Tables II, III and IV. Second, PARAFAC iden- 78 tifies the amplitude and phase providing channels individually rather than identifying the actual channel pairs. This results in P × Q possible channel pairs, where P is the number of amplitude providing channels and Q is the number of phase providing channels. In order to determine the actual coupled channel pairs, one needs to do significance testing which is computationally expen- sive. Estimation of coupled channel pairs from the indices of the sparse matrix in the HoRPCA method overcomes these limitations. Based on the F-measure values reported for various simula- tions, HoRPCA based method outperforms PARAFAC and averaging based method and is more robust to varying signal conditions such as noise level, subject variability and number of cou- plings. HoRPCA is also more robust to noise compared to gedCFC indicating the advantage of tensor representation with respect to matrix factorization based methods (Table 3.1). In the analysis of EEG data, the proposed approach detected event-related theta-gamma PAC during error and correct responses with higher PAC for error response. This finding is consistent with prior studies where theta-gamma coupling was reported during visual tasks like working memory processing and serial memory recall [52]. Theta-gamma PAC was also reported in an error processing MEG study [111] and in error/correct trial learning [21]. A possible explanation for the presence of higher theta-gamma PAC during error response is that error trials may reflect a miscoding of information, which leads to a large-scale functional integration across theta-gamma frequency bands to improve the performance after an error response [21]. Cross-frequency networks extracted through HoRPCA are sparse networks as the background couplings are discarded (Figure 3.3, Figure 3.6). Analysis of the low-rank parts of the tensors for both error and correct responses showed that there was no significant difference between the two response types based on solely the low-rank tensor (Figure 3.3). This observation justifies our hypothesis that the low-rank part of the tensor captures background PAC activity while the sparse part corresponds to response-evoked PAC. Observation of the detected PAC networks shows that 79 the networks corresponding to the correct response were concentrated between frontal theta phase providing channels and parietal gamma amplitude channels. On the other hand, networks for error response were concentrated between frontal-central theta phase activity and parietal gamma activ- ity. Comparison of ERN and CRN networks through statistical significant testing revealed that the difference was centered around the medial frontal cortex, which has been previously reported to be active during error processing [111]. Prior studies also hypothesized that error-related negativity initiates the medial frontal based top-down control mechanisms to improve the performance after an error response [96, 97]. Thus, the PAC networks extracted by the proposed HoRPCA approach are consistent with previous literature reflecting higher theta-gamma coupling in the medial frontal cortex and relating this with error-related negativity. The high classification accuracy achieved by HoRPCA based multivariate PAC also indicates that multivariate PAC can be successfully em- ployed to characterize cross-frequency connectivity networks associated with error and correct responses. Results from interrelation study revealed estradiol was associated with an enhancement of error-related theta-gamma coupling both within females across the menstrual cycle and between females with varying levels of average circulating estradiol. Findings suggest engagement of a frontoparietal network in error monitoring that is modulated by the endocrine system in females. Together with prior work, the present results further show the importance of inter-areal PAC study that revealed the role of midfrontal theta in modulating higher frequency power bursts across cortex that are involved in cognitive control-related brain functions. 3.7 Conclusion In this study, we proposed a HoRPCA based multivariate PAC framework to capture the dynam- ics of cross-frequency PAC across various brain regions and frequencies. The empirical results 80 obtained from the analyses of simulated and EEG data supported our hypothesis that the sparse tensor extracted through HoRPCA can capture the spectral and spatial dynamics of event-related multivariate PAC network while discarding the background perturbations as part of the low-rank tensor. With these unique properties, the proposed HoRPCA based multivariate approach can lead to obtaining a complete understanding of connectivity across frequency bands and brain regions. 81 Chapter 4 Assessment of Dynamic Phase Amplitude Coupling Using Matching Pursuit 4.1 Introduction Neural oscillations in the brain play a fundamental role in cognitive processing by demonstrat- ing coordinated interactions among different frequencies [132, 133, 14, 9]. The interactions across multiple frequencies are referred to as Cross Frequency Coupling (CFC) [14, 13, 19]. CFC exhibits itself in many forms including phase to phase [9], amplitude to amplitude [25], phase to frequency [13] and phase to amplitude coupling [14, 52]. Phase amplitude coupling (PAC) has gained a lot of attention as a CFC measure. PAC refers to the modulation of the amplitude of the higher frequency oscillation by the phase of the lower frequency oscillation. Numerous empirical findings and theo- retical predictions link PAC with neural coding and information transfer within the complex neural networks [14, 19, 52, 65]. This phenomenon was observed in various neurophysiological signals collected from various species and within multiple brain regions during multiple tasks [52]. Recent research also relates disruption of PAC with pathological brain states like autism spectrum disor- der [134, 135], schizophrenia [136], and Parkinson’s Disease [137]. Therefore, an unbiased and accurate estimation of PAC is necessary for investigating the adaptive and pathological elements of cross-frequency neural interactions. 82 The conventional method for quantifying PAC is based on identifying the amplitude of the high frequency and the phase of the low frequency oscillation [14, 52]. As such, in the first step, the sig- nals are bandpass filtered to extract the low and high frequency oscillations [14, 52]. An analytical signal is then obtained for each frequency using the Hilbert transform, from which instantaneous phase and amplitude components are obtained to compute PAC [52]. PAC is then quantified using these instantaneous phase and amplitude components with metrics such as the mean vector length [14], modulation index [52] and the phase-locking value [138]. Although this Hilbert transform based conventional PAC has been successfully employed, it has several drawbacks. First, this measure is sensitive to the choice of the bandpass filter parame- ters [139]. If the selected bandwidth is too broad or narrow, it is prone to result in false positives and negatives, respectively [140]. Second, this measure is stationary in nature, focusing on phase amplitude modulations within a certain time period or over arbitrary sliding short time windows. As brain interactions are highly dynamic [141], it is necessary to characterize the temporal dy- namics of neural networks with millisecond of accuracy [142]. As such, a reliable estimation of PAC requires not only identifying the existence of coupled frequencies but also determining the variation of PAC across time. A few methods have been proposed to capture the dynamics of PAC. The most common approach is to compute PAC using an arbitrary time window [52, 143]. More recently, sliding time window based PAC computation was introduced in [80]. However, this method requires the user to choose the window length and amount of overlap, which substan- tially affects the quantified PAC. An event-related phase-amplitude coupling (ERPAC) measure was introduced by Voytek et al. to evaluate the task-related changes in PAC across events [144]. However, this method cannot provide a trial-by-trial description of PAC variability and needs a collection of experimental events that are assumed to generate equivalent signal dynamics. In this work, we propose a novel dynamic PAC computation approach that addresses these lim- 83 itations. This procedure is based on matching pursuit (MP). First, we decompose the signal into time and frequency localized atoms using an enriched dictionary [145] which includes the Gabor functions along with the arctangent function. The enriched dictionary of atoms allows flexible and sparse representation of the signal resulting in a more accurate description of signal dynamics. Using the properties of these time and frequency localized MP atoms, we estimate the coupled low-high frequency pairs by employing k-means clustering. We quantify PAC by extracting the instantaneous phase and amplitude components from these atoms, which allow us to capture the temporal variation of PAC using a data-driven approach. Moreover, unlike the conventional Hilbert transform based methods, this approach does not require any bandpass filtering. We evaluated the performance of the proposed measure on both simulated and event-related EEG data. The simu- lated signals contain PAC at different low-high frequency combinations along with noise and vari- ation of coupling duration in the signal. We also compare the performance of MP based dynamic PAC with the conventional sliding window based method. 4.2 Methods 4.2.1 Phase Amplitude Coupling (PAC) For a given signal x(t), PAC between the high frequency ( fa ) and low frequency ( f p ) oscillations is defined as the modulation of the high frequency amplitude, A fa (t), by the phase φ f p (t) of the low frequency oscillation. The computation of PAC is performed as follows: 1. x(t) is bandpass filtered at the two frequencies of interest (high frequency fa and low fre- quency f p ). These filtered signals are denoted as x fa (t) and x f p (t), respectively. 2. The amplitude component A fa (t) at fa is extracted from the Hilbert transform of x fa (t). The 84 Hilbert transform (H(t)) ˜ of the filtered high frequency signal is computed as [146]: Z+∞ 1 x fa (τ) H˜fa (t) = P.V. dτ, (4.1) π t −τ −∞ where P.V. indicates that the integral is taken in the sense of Cauchy principal value. The amplitude component A fa (t) is then extracted as: A fa (t) = |H˜fa (t)|. 3. The phase component φ f p (t) at f p is extracted from the Hilbert transform H˜f p (t) of x f p (t) as:   H˜f p (t) φ f p (t) = arg |H˜ (t)| . fp 4. After extracting the amplitude and phase components, the phase-amplitude coupling between these two components is computed using the mean vector length (MVL) measure which quantifies the coupling between f p and fa by taking the length of the average vector [12] as follows: 1 T MV L( fa , f p ) = ∑ A fa (t)e jφ f p (t) , (4.2) T t=1 where T is the number of time points. 4.2.2 Sliding Window based Dynamic PAC In this study, we compare our proposed approach with previously proposed sliding window based dynamic PAC approach. This approach computes PAC over time windows that slide along the length of the given signal x(t). Parameters like the length of sliding time window and overlap between these windows are defined by the user. The amplitude and phase components are extracted using Hilbert transform within the given time window and PAC values are computed using (4.2) to 85 assess dynamic changes in PAC. The sliding window length is chosen such that it is long enough to cover at least one full cycle of the lowest frequency of the low frequency band of interest and is moved along the data time length with a 50% of overlap [80]. As suggested in [80], to improve the frequency resolution and minimize the filter edge effects, the filtering was performed on the full length signal before applying the sliding windowing. 4.2.3 Matching Pursuit Matching pursuit [147] is a commonly used greedy algorithm which aims to obtain a sparse linear representation of a signal, x(t), in terms of functions (also referred to as atoms), gi for i ∈ {1, 2, . . . , P}, from an overcomplete dictionary, D. MP provides high-resolution signal repre- sentation and allows for parametric description of both periodic and transient signal features. As finding an optimal approximation of a signal by selecting functions from a dictionary is NP-hard, a suboptimal iterative search algorithm [147] is used as follows: 1. Define the 0th order residual as R0 x = x, D̃ = D and set k = 0. 2. For the kth order residual, Rk x, the best atom is selected such that the inner product between the atom and the residual is maximized as: gk = argmaxgi ∈D̃ |hRk x, gi i|. 3. Compute the (k + 1)th order residue as: Rk+1 x = Rk x − hRk x, gk igk . 4. Set k = k + 1, D̃ = D̃ \ gk , and go back to step 2 until a pre-determined number of atoms is achieved or the normalized reconstruction error is below a pre-determined threshold such as k k2  Rk+1 x  kxk < 0.05 . 2 86 After M iterations, the following linear representation is obtained: M x= ∑ hRk x, gk igk + RM+1x. (4.3) k=1 This procedure converges to x in the limit, that is, x = ∑∞ k k=1 hR x, gk igk , and preserves signal energy. 4.2.4 Matching Pursuit based Dynamic PAC (MP-dPAC) In this chapter, we propose a MP based dynamic PAC measure denoted as MP-dPAC. Using an enriched dictionary proposed in [145] which includes functions with varying degrees of asymmetry – from symmetric Gabor atoms to highly asymmetric waveforms. This dictionary allows us to capture the presence of asymmetric waveforms in signals. The functions of this dictionary are constructed in two parts: the ascending part is based on a Gabor function, and the descending part on an arctangent function. The elements of this dictionary can be described as: −(t−µk )2 ! 2σk gk (Nk , αk , β , µk , σk , fk )(t) = Nk exp atan(β (t−µk ))+ π2 e j2π fk t , (4.4) 1 + αk (t − µk ) π where Nk is the normalization constant such that ||gk ||= 1, σk is the scale or time span parameter, µk is the time center, fk is the modulation frequency, αk = 1/2σk 2 is the slope of the descending function and β is a constant that is high enough to capture the whole range of atan argument (here we used the default β = 1016 ). For a signal with T = 2L samples, where L ∈ Z, the scale σk is chosen from the dyadic sequence in terms of octave j (integer) such that σk = 2 j and parameters µk and fk are sampled dyadically within the interval. To compute the MP based dynamic PAC of a signal, x(t), first the signal is decomposed into M 87 atoms using the implementation of MP algorithm in EEGLAB Matlab Toolbox [148]. The signal was decomposed until normalized reconstruction error is less then 5% or M = 20 whichever came first. The maximum number of atoms was set to 20 based on the prior work in EEG literature [149, 150] that indicates that Gabor atoms (and their frequency modulated versions) are very effective at compressing EEG signals with a few number of atoms. The time, frequency centers and the scale of the atoms are then used to reveal the dynamics of the underlying signal [145]. The steps for computing MP based dPAC are as follows: 1. Identifying Pairs of Coupled Atoms: After extracting the MP atoms, our objective is to identify the atoms that are coupled to each other during a certain time span. With this objective, for each atom gi , we obtain the start time (tstart ), stop time (tstop ), time span (tspan ) and frequency ( fi ). To identify pairs of atoms that are operating within the same time window, we perform k- means clustering with the (tstart ), (tstop ) and (tspan ) as the features describing each atom. k-means clustering with squared Euclidean distance is performed following Lloyd’s algorithm [151], to assign the atoms into exactly one of c ∈ {1, 2, ..,C} clusters defined by centroids. The optimal number of clusters is determined using the elbow method [152]. The atoms assigned to the same clusters are considered to be potentially coupled as they are active during the same time span. 2. Quantify PAC between the Coupled Atoms: Once the atoms that are in the same cluster are detected, we need to verify if there is any coupling between these atoms. First, we identify the low frequency ( f p ) and high frequency ( fa ) atoms in the same cluster c. An atom was selected as low frequency if the center frequency of the atom was within one of these low-frequency bands: delta (1-3 Hz), theta (4-7 Hz), alpha (8-12 Hz), and beta (13-30 Hz), and was selected as high frequency if the atom frequency was within gamma (31-100 Hz) band. After selecting the high and low frequency atom pair within a cluster, we extract the amplitude component from the high frequency atom and phase component from the low frequency atom. As these atoms are time and 88 frequency localized, there is no filtering involved unlike the PAC procedure described in Section 4.2.1. For a pair of atoms [gcflp , gcfma ] located in cluster c, with l and m indicating the atom indices, f pcl is the center frequency of gcflp and facm is the center frequency of gcfma . Hilbert transform Hgcm (t) fa of the high frequency atom gcfma is computed following (4.1) and amplitude component Acfma (t) is extracted as: Acfma (t) = |Hgcm |. (4.5) fa Next, the phase component at f pcl , is extracted by computing the Hilbert transform Hgcl (t) of fp the low frequency atom gcflp following (4.1), and the phase component φ cf pl (t) is extracted as:  H cl (t)  g fp φ cf pl (t) = arg . (4.6) |Hgcl (t)| fp After detecting the amplitude and phase components, the phase-amplitude coupling between these two components is computed using the MVL measure as follows: c 1 cm jφ f pl (t) MV L( f pcl , facm )(t) = ∑ fa A (t)e , (4.7) |Tc | t∈T c where c is the cluster of the selected atom pairs and |Tc | is the number of time points in the time window Tc = Tcl ∩ Tcm , Tcl = [tstart (gcflp ),tstop (gcflp )] and Tcm = [tstart (gcfma ),tstop (gcfma )]. 3. Permutation Testing: To ensure the validity of the computed dynamic MP-dPAC, statistical significance testing was performed by generating surrogate datasets. Following the guidelines in Aru et al. [114], we generated surrogate datasets using a block swapping approach. A permuted 89 Figure 4.1: Illustration of the proposed matching pursuit(MP) based dynamic PAC method. Steps 1 to 9 highlight the methodologies followed for computing the MP based dynamic PAC. signal was generated by splitting the amplitude envelope of the high frequency oscillation, Acfma (t) at a random time point and then swapping the two resulting time series. MVL values are then computed by associating the permuted amplitude Acfma (t) of high frequency oscillation with the original phase φ cf pl (t) of the low frequency oscillation. This whole procedure was performed 500 times to generate a sequence of surrogate MVLs observed under the null hypothesis that phase- amplitude coupling is due to chance. We compute 95% confidence intervals for the surrogate distribution. The observed MP-dPAC value was deemed to be significant only when it is larger than 95% of the surrogate MP-dPAC values. The observed coupling value is normalized based on MV Lobserved −µMV Lsurr the distribution of the surrogate MVLs by computing the z value as MV Lz = σMV Lsurr , where µMV Lsurr is the mean and σMV Lsurr is the standard distribution of the surrogate distribution. 90 The methodological steps detailed above are summarized in Figure 4.1. We repeat the whole procedure for all the atoms in the same cluster and then for all clusters. 4.2.5 Synthesized Data Synthesized data with dynamic coupling were generated by varying the modulation of the high frequency amplitude signal by the low frequency phase signal. A three second signal was generated to introduce dynamic coupling. For the first one second, the phase of slow oscillation at f p 1 = 15Hz was coupled to the amplitude of a faster oscillation at fa 1 = 50Hz. For the last two seconds, the first coupling mode was terminated and two other modes were generated simultaneously with f p 2 = 10Hz coupled with fa 2 = 60Hz, and f p 3 = 5Hz coupled with fa 3 = 40Hz, respectively. The modulated signal was generated as:   0 ≤ t < 1,  x f p (t) + x fa (t),  1 1 x(t) = , (4.8)   x f p (t) + x fa (t) + x f p (t) + x fa (t), 1 ≤ t < 3  2 2 3 3 where the low frequency phase signal was generated as x f p (t) = K f p sin(2π f pt), with K f p = 1 and f p ∈ { f p1 , f p2 , f p3 }. The high frequency amplitude signal was generated as, x fa (t) = A fa (t) sin(2π fat), jφ f p (t) ) where, A fa (t) = K fa 1+Re(e2 , and K fa ∈ [0, 1] determines the coupling strength, φ f p (t) is the phase of the low frequency phase providing signal and fa ∈ { fa1 , fa2 , fa3 }. The coupling strength was kept constant and identical in all modes. Following [80], additive noise was added to the signal which was generated by combining white Gaussian noise with (1/ f ) noise. The power law sam- ples simulate the background brain activity, and the white Gaussian samples simulate measurement noise. The variance of the white Gaussian noise was set to half of the variance of the power law  samples. Signal to Noise Ratio (SNR), is defined as the ratio of signal power Psignal to the noise 91   Psignal power (Pnoise ) and is expressed in decibels as SNRdB = 10 log10 Pnoise . 4.2.6 EEG Data An EEG dataset collected during a cognitive control-related error processing study was used to assess the proposed dynamic PAC measure. The experiment performed was a letter version of the speeded-reaction Flanker task [128]. The study was designed following the guidelines and experimental procedure approved by the Institutional Review Board (IRB) of the Michigan State University. Informed consent for experimentation was obtained from each participant before data collection. For each trial of the Flanker task, a series of five letters, either congruent (e.g., AAAAA) or incongruent (e.g., AABAA), were displayed in front of each participant. The participants iden- tified the center letter using a mouse regarding the Flanker letters. Each trial started with a flanking stimulus (e.g., AA AA) of 35ms. Then, the target stimuli were inserted in the middle of the flanker letters (e.g., AAAAA/AABAA) and displayed for 100 ms, which resulted in 135ms presentation time. An inter-trial interval of variable duration (1200 to 1700 ms) was then inserted. Continuous EEG data were recorded during this task using the ActiveTwo system (BioSemi, Amsterdam, The Netherlands). 64 Ag–AgCl channels were placed following the international 10/20 system. The sampling frequency of the data was 512 Hz. The trials with artifacts were re- moved. The volume conduction was minimized using the Current Source Density (CSD) Toolbox [90]. The pre-processed data were used for the analysis. The study focused on trials corresponding to Error-related negativity (ERN) occurring after an error response and the Correct-related nega- tivity (CRN) occurring after a correct response. The length of each trial was one second. A total of 34 participants were considered for the ongoing study. As the total number of trials corresponding to correct/error response varies for different participants, the number of trials considered for both 92 responses was kept the same (20 trials) for a fair comparison. 4.3 Results 4.3.1 Synthesized Data Results The performance of the proposed MP based dynamic PAC was evaluated in terms of both the accuracy of detecting the coupled frequency pairs and estimation of the coupling duration. The input signal was generated using the model described in Section 4.2.5. The clustering of the atoms using the (tstart ), (tstop ) and (tspan ) parameters as feature resulted into two clusters as shown in Figure 4.2. The coupling strength of the respective atom pairs was computed using (4.7). Using these MVL values, comodulograms were generated to exhibit the cross-frequency coupling between different oscillations, with low frequency values [2, 20] Hz along the x-axis, and high frequency values [30, 80] Hz along the y-axis. Figure 4.3 (a) depicts the comodulogram constructed by the proposed MP-dPAC as a typical, non-time-resolved representation of coupled frequencies. To further illustrate MP-dPAC’s ability to track the time-varying PAC, time varying coupling strength over low and high frequency ranges are shown in Figures 4.3 (b) and (c), respec- tively. We also show similar plots constructed using SW based approach in Figures 4.3 (d), (e) and (f). Comparing the comodulograms from both methods (Figures 4.3 (a) and (d)), it can be seen that the proposed method correctly detects the phase (5 Hz, 10 Hz and 15 Hz) and amplitude (40 Hz, 60 Hz and 50 Hz) providing frequencies with higher resolution. Unlike the SW based method, the proposed method does not require bandpass filtering of the signals resulting in high resolu- tion comodulograms. The coupling strength vs frequency maps reveal the time course of all three couplings in the data for both approaches. Comparing these figures, we can see that MP-dPAC can capture the dynamic nature of PAC more precisely. It can also identify the frequencies more 93 accurately whereas with SW there are some additional frequency components. Figure 4.2: Clustering performance of the extracted MP atoms from synthesized data using the (tstart ), (tstop ) and (tspan ) parameters. To evaluate the effect of coupling duration on the proposed MP-dPAC, we performed a system- atic examination of the sensitivity of PAC estimates to coupling duration. Signals were simulated with a theta-gamma (5 Hz coupled with 60 Hz) coupling where the length of the coupling duration was varied. The length of the signal was fixed at two seconds with a sampling frequency of Fs = 1000 Hz. The coupling duration was varied in terms of T s = 1/Fs from 0.2/T s = 0.2 second to 1.6/T s = 1.6 second. SNR was fixed at 6dB. We compare the performance of MP-dPAC with the SW based dynamic PAC method [80]. Figure 4.4 summarizes the performance of the tested meth- ods in terms of PAC value (mean ± std) for 100 simulations for each time duration. As depicted in Figure 4.4, the proposed MP-dPAC measure can detect the dynamic coupling more accurately compared to sliding time window based method for time durations as low as 0.2/T s (200 millisec- 94 Figure 4.3: Comparison of Matching Pursuit (MP) based dynamic PAC approach with sliding window (SW) based dynamic PAC approach for synthesized data. (a) Signal with dynamic coupling; (b) MP-dPAC comodulogram ( f p vs. fa map); (c) MP-dPAC coupling strength map over time for low frequency ( f p vs. time map); (d) MP-dPAC coupling strength map over time for high frequency ( fa vs. time map); (e) SW-dPAC comodulogram ( f p vs. fa map); (f) SW based PAC coupling strength map over time for low frequency ( f p vs. time map); (g) SW based PAC coupling strength map over time for high frequency ( fa vs. time map). The signal to noise ratio was fixed at 6dB. 95 ond). We also compare the methods by quantifying the slope of PAC versus coupling duration curves. The average slope computed across different coupling durations is 1.147 for MP-dPAC, compared to 2.932 for SW-based dPAC. We also tested the effect of varying noise level on the dynamic PAC estimation performance. Synthesized signals were generated with theta-gamma coupling (5 Hz coupled with 60 Hz).The length of the signal was one second and the coupling duration was 800ms. SNR values were varied from -10dB to 10dB. Comparison of SW to MP based dynamic PAC was performed and is reported in terms of (mean ± std) for 100 simulations for each SNR value as shown in Figure 4.5. The average slope for MP-dPAC with respect to SNR is 0.1231 compared to 0.4761 for SW based dPAC. 20 std) 18 16 Dynamic PAC (mean 14 12 10 8 MP-dPAC SW-dPAC 6 0.2/Ts 0.4/Ts 0.6/Ts 0.8/Ts 1.0/Ts 1.2/Ts 1.4/Ts 1.6/Ts Duration of coupling (s) Figure 4.4: Performance comparison of dynamic PAC values for MP-dPAC and sliding win- dow based dPAC for different coupling durations. Here, T s = 1/Fs and Fs = 1000. SNR was fixed at 6dB. These results demonstrate that the proposed measure is more robust to noise and variation of 96 20 std) 15 Dynamic PAC (mean 10 5 MP-dPAC SW-dPAC 0 -10 -5 0 5 10 Signal to Noise Ratio (dB) Figure 4.5: Performance comparison of dynamic PAC values for MP-dPAC and sliding win- dow based PAC for different noise levels. coupling duration compared to sliding window based method. Moreover, the standard deviations reported in Figures 4.4 and 4.5 show that the variability of the proposed MP-dPAC method is lower than the SW based method across multiple realizations of the experiments. This also highlights the stability of the proposed method. 4.3.2 EEG Data Results MP-dPAC was computed for all 64 channels for error and correct response for all 34 subjects. First, we compute the number of atoms per trial, channels, clusters along with the number of significantly coupled atoms per cluster for both response types averaged across all channels and subjects as listed in Table 4.1. It should be noted that each channel has 20 trials and the atoms extracted from all trials are combined before clustering to compute PAC across the trials. The number of atoms per trial averaged across subjects and channels for each response is lower than 97 10 with a small standard deviation which shows that the representation of the EEG signal provided by this enriched dictionary is sparse and does not vary much across different trials and subjects. Similarly, the number of clusters is between 3 ∼ 4 indicating a small number of time windows of interest. While the number of atoms per cluster is high, the number of significantly coupled atoms per cluster is around 13 ∼ 15 range indicating the different coupling modes within that time window. For the analysis of EEG data, the low frequency ( f p ) and high frequency ( fa ) ranges of interest were defined as [2, 12] and [40, 80] Hz, respectively. We focused on a time range of -200ms before the response to 400ms after the response. We refer to -200ms to 0ms as pre-response, 0- 200ms as response and 200ms-400ms as post-response time windows. Figure 4.6 shows the time varying coupling strength averaged across subjects over low and high frequency ranges for error and correct responses. Figure 4.6(a) and (b) show the time varying coupling strength for error response whereas Figure 4.6(c) and (d) show the time varying coupling strength averaged across the subjects over low and high frequency ranges for correct response. From these figures, it can be seen that the MP-dPAC averaged across subjects was higher for ERN compared to CRN. For ERN, highest PAC values were concentrated around FCz whereas for CRN the highest PAC values were concentrated around F6. EEG data analysis also revealed that MP-dPAC for theta-gamma and alpha-gamma bands is more prominent for ERN compared to CRN and these differences are observed mainly during response and post-response time windows. As FCz provided the highest MP-dPAC value for error response, the coupling strength at FCz is depicted both as comodulogram ( f p vs. fa ), and coupling strength over time map in terms of f p vs. time and fa vs. time, for both error and correct responses in Figure 4.7. Similarly, as F6 has the highest MP-dPAC for correct response, we show the coupling strength at F6 both as comodulogram and coupling strength over time map for both error and correct responses in Figure 4.8. We also 98 Table 4.1: The number of MP atoms per trial, electrodes, clusters along with the significantly coupled atoms per clusters for error and correct response types averaged across all electrodes and subjects. Parameter Error response Correct response Number of atoms per trial 8.815 ± 0.451 9.517 ± 0.567 Number of atoms per electrodes 194.580 ± 24.486 201.328 ± 28.574 Number of clusters 3.639 ± 0.766 3.830 ± 1.384 Number of atoms per clusters 44.647 ± 6.991 40.385 ± 6.957 Number of significantly coupled atoms per clusters 15.579 ± 3.785 13.840 ± 3.184 show the comodulogram, coupling strength over time map in terms of f p vs. time and fa vs. time for the SW based approach in the third and fourth columns of Figures 4.7 and 4.8. From these figures, it can be observed that comodulograms for these channels indicate a strong theta-(low) gamma coupling for ERN along with some alpha-(high) gamma coupling for both ERN and CRN. Coupling strength over time map tracks the dynamics of PAC in these channels. Figure 4.7 shows that the FCz channel exhibited highest average theta-gamma coupling (4-7 Hz coupled with 42∼48 Hz) for ERN and this coupling is concentrated mainly around the response time win- dow (17.74 ± 11.89ms to 192.63 ± 23.46ms). The alpha-gamma coupling (8-12 Hz coupled with 72∼80 Hz) for ERN was concentrated around post response time window (153.61 ± 33.12ms to 278.25 ± 43.69ms). Thus, MP-dPAC reveals that the theta-gamma coupling occurs at the start of ERN response, whereas the alpha-gamma coupling spreads beyond the response time and lasts longer than the theta-gamma coupling for FCz. As shown in Figure 4.8, for F6, alpha-gamma coupling is concentrated around post-response time window for CRN while it spreads from pre to post-response time window for ERN. Comparing the PAC detected from MP and SW methods, we can see that both methods detected similar sets of coupled low-high frequency pairs. However, comparison between the coupling strength over time maps shows that the SW based method fails to capture the coupling dynamics. The observed MP-dPAC values averaged across all channels and subjects for error and cor- 99 Figure 4.6: Coupling strength over time MP-dPAC map for error (first row) and correct (second row) responses. (a-b) MP-dPAC coupling strength over time map during error response for (a) low frequency ( f p vs. time map) and (b) high frequency ( fa vs. time map); (c-d) MP-dPAC coupling strength map over time map during correct response for (c) low frequency ( f p vs. time map) and (d) high frequency ( fa vs. time map). rect response types for three different time windows (pre-response, response and post response) are listed in Table 4.2. This table depicts that the average theta-gamma coupling was highest for the ERN time window whereas the average alpha-gamma coupling was highest for post-response time window. We use the dynamic PAC computed from different time windows for all subjects 100 Figure 4.7: Dynamic PAC for EEG data during error and correct response for FCz channel in terms of comodulogram and coupling strength over time map. (a, b, c) Dynamic PAC during error and (d, e, f) dynamic PAC during correct response for FCz channel in terms of comodulo- grams (a-d), and coupling strength over time map for low frequency (b-e) and high frequency (c-f) computed using matching pursuit (first and second column) and sliding window (third and fourth column) based approaches. to determine the brain regions which exhibited the most differentiation between error and correct responses. A Wilcoxon Signed Rank Sum Test with Bonferroni correction is conducted to deter- mine the statistical significance of the difference between ERN and CRN for both theta-gamma and alpha-gamma coupling across subjects for different time windows. Figure 4.9 shows the p-value topo plots that depict the channels with significant differences between ERN and CRN corre- sponding to theta-gamma (first row) and alpha-gamma (second row) couplings. For theta-gamma coupling, the significant difference between the two response types was highest during response time window (0ms-200ms) with p = 0.0252 and was concentrated around frontal-central and pari- etal brain regions. For alpha-gamma PAC, even though there are some significant PAC differences during response time window in the same regions, the level of significance is p = 0.0412 and the regions that show significant differences are localized around fewer channels. The alpha-gamma 101 Figure 4.8: Dynamic PAC for EEG data during error and correct response for F6 channel in terms of comodulogram and coupling strength over time map. (a, b, c) Dynamic PAC during error and (d, e, f) dynamic PAC during correct response for F6 channel in terms of comodulo- grams (a-d), and coupling strength over time map for low frequency (b-e) and high frequency (c-f) computed using matching pursuit (first and second column) and sliding window (third and fourth column) based approaches. coupling was found to be more significantly different (p = 0.0316) for error vs. correct during post-response time window (200ms-400ms) and this difference was concentrated around frontal and occipital brain regions. Table 4.2: The MP-dPAC averaged across all channels and subjects for error and correct response types for different time windows. Time Theta-Gamma PAC Alpha-Gamma PAC window Error Correct Error Correct Pre-response (-200-0 ms) 3.75± 1.05 3.56± 1.21 3.75± 1.14 3.34± 1.18 Response (0-200 ms) 5.41 ± 1.02 3.84 ± 1.23 3.93 ± 1.08 3.67 ± 1.17 Post-response (200-400 ms) 4.26 ± 1.04 3.79 ± 1.12 5.07 ± 1.06 3.86 ± 1.15 102 Figure 4.9: p-value topo-plots showing significant difference between error and correct re- sponses for theta-gamma (first row) and alpha-gamma (second row) coupling. Column 1 (a and d) shows the topoplots for pre-response (-200ms to 0ms) time window, column 2 (b and e) shows the topoplots for response (0ms to 200ms) time window and column 3 (c and f) shows the topoplots for post-response (200ms to 400ms) time window. The color scale from white to black indicates the level of significance (Wilcoxon Signed Rank Sum Test with p < 0.05). 4.4 Discussion Human brain routes information through a multitude of oscillatory frequencies [133], and un- derstanding how the interactions occur across various frequency bands is a fundamental research question [14, 11]. PAC has been reported as one of the most prominent cross-frequency coupling mechanisms that controls the neural integration across populations of neurons [14]. Many methods exist to quantify PAC, each with specific merits, but most approaches assume PAC is stationary in nature. As brain activity is highly dynamic [141], ideal PAC measure should be able to capture the dynamic nature of PAC. In this article, we have proposed a time-varying PAC measure for 103 quantifying dynamic neural coupling in the brain. The proposed measure is data-driven, based on matching pursuit algorithm. MP was implemented using an enriched dictionary of functions of various degrees of asymmetry – from symmetric Gabor atoms to highly asymmetric waveforms [145]. This approach provides a flexible and sparse representation of the time series, allows for the correct determination of the time-span of the components and removes the “pre-echo” or “energy leakage” effect which may cause the appearance of energy even before the start of the signal [145]. It is crucial to remove the “energy leakage” effect to obtain all the atoms that describe the signal properly rather than incorrectly extracting one long-lasting atom from the signal. The effectiveness of the proposed MP-dPAC method was verified for both simulated and neu- rophysiological data recorded from human participants. The results from synthesized data (Figure 4.3) showed that the proposed MP-dPAC method detects all the coupled frequencies and the time variation of the coupling correctly with high time and frequency resolution. This is because this method estimates both the high frequency amplitude component and low frequency phase com- ponent from the time and frequency localized atoms. Using the properties of these atoms, the proposed method can extract high resolution estimates of both the envelope and the phase of the oscillations resulting in high resolution comodulogram. Traditional Hilbert transform based meth- ods rely on bandpass filtering and it is well established in the literature that the filter properties substantially hamper the coupling estimation [114, 153, 154]. Since MP-dPAC does not require any filtering, it alleviates the limitations associated with this filtering step. The proposed MP-dPAC measure is also shown to be more effective at capturing PAC with a short duration (Figure 4.4) and provide improved robustness under noise (Figure 4.5) compared to the SW based approach. This is because PAC value computed using the SW based approach is an average index computed across a pre-determined arbitrary time window and as such, does not adapt to the signal dynamics. A window length that is selected too broad or narrow may result in an incorrect measure of PAC. For 104 MP-dPAC, on the other hand, the MP decomposition parameterizes the time-frequency surfaces and by computing PAC using these time and frequency localized MP atoms, we overcome the limitations of window length selection. MP decomposition also offers a succinct representation of the data and thus acts as an effective denoising method [149]. Therefore, the proposed data-driven MP-dPAC approach offers improved performance at capturing the dynamic nature of PAC over the SW based approach. In the EEG data analysis, the proposed approach detected theta-gamma PAC during the re- sponse window for error response and alpha-gamma PAC during the response window and post- response window for both error and correct responses. The coupling strength over time maps depict that the MP-dPAC analysis disentangles the theta-gamma and alpha-gamma coupling and indicates that the theta-gamma coupling is shorter in duration and is more associated with ERN. In contrast, though overlapping in time, the alpha-gamma coupling is longer in duration and is not di- rectly related to ERN. This information could not be obtained through traditional comodulograms. Higher PAC values were observed for error response compared to the correct response. This find- ing is consistent with prior studies reporting higher theta-gamma PAC in an error processing MEG study [111] and error/correct trial-learning study [21]. An increase in PAC was also associated with working memory processing [110] and serial memory recall visual tasks [52]. The increase in PAC for error response fortifies the idea that the error trials may reflect a miscoding of information, which leads to a large-scale functional integration across different frequency bands to improve the performance after an error response [21]. Prior studies also reported an increase in the number of regions engaged together in performance monitoring [155]. For ERN, the computed MP-dPAC values are concentrated around the central channels and most prominent for FCz. This region was shown to be active during negative feedback and error-related studies [111, 156, 157] and it was suggested that the error response initiates the medial frontal based top-down control mechanisms 105 to improve the performance [96, 97]. Moreover, the p-value topoplots (Figure 4.9) summarize the dynamic evolution of event-related PAC. From our extensive experiments, we have observed that the proposed method is robust to the different choice of the design parameters (coupling duration and noise). In particular, as MP con- verges very quickly for EEG data, the number of atoms is usually low, so our atom selection criterion does not affect the results. The choice and the size of the dictionary are important. As long as the selected dictionary is an overcomplete collection of time-frequency localized atoms, the results will not vary much. Moreover, the choice of the clustering method is not key to the accuracy of the algorithm as the current implementation relies on a simple k-means clustering al- gorithm. Finally, similar to other univariate measures of PAC, the proposed measure is limited in the sense that it can only quantify PAC within a brain region rather than quantifying PAC across different brain regions. As most physiological data offers multi-channel recordings, in future re- search, the proposed MP-dPAC framework can be extended to multivariate settings for all possible frequency band combinations to capture the dynamics of PAC across time, frequency, and spatially distributed neuronal locations. 106 Chapter 5 Conclusions and Future Work 5.1 Conclusions In this thesis, we proposed several methods to detect and quantify the phase amplitude coupling. The proposed methods have significant contributions to the field of cross-frequency coupling in neuronal oscillations. In order to better quantify PAC in complex brain networks, we started our work by proposing a new PAC measure and then generalized it to analyze inter-areal PAC between different brain regions. After that, we presented a dynamic PAC measure to track the variation of PAC with time. The main focus of this thesis is to capture the spatial, spectral, and temporal dynamics of cross-frequency phase-amplitude coupling in the brain. In Chapter 2, we introduced a novel time-frequency based PAC method for estimating cross- frequency phase-amplitude coupling and illustrate that the proposed method offers higher accuracy and robustness compared to existing methods through simulations. The proposed time-frequency based PAC method provides several advantages over existing Hilbert and Wavelet-MVL methods. The proposed method estimates both the amplitude of the high frequency signal component and the phase of the low frequency component directly from a high resolution complex time-frequency distribution. Using properties of Cohen’s class of distributions, the proposed method can extract high resolution estimates of both the envelope and the phase of the oscillations resulting in high resolution comodulograms. Unlike Hilbert transform based methods, the proposed method does 107 not require bandpass filtering of the signals, and thus overcomes the problems of misidentification of the frequencies due to the selection of the bandwidth and the transition band of the bandpass filter. While analytic wavelet transform based PAC methods address some of the shortcomings of Hilbert transform based PAC estimates, they are highly dependent on the selected wavelet parame- ters and the user provided range of frequencies. Unlike wavelet based PAC measures, the proposed approach does not depend on any design parameters and offers uniform frequency resolution. Thus the estimated PAC comodulograms do not change with the range of frequencies considered. The proposed algorithm was evaluated on multiple synthetic and real neurophysiological signals. The results show that the proposed PAC measure is more robust to varying signal conditions such as noise level, signal length, coupling strength and the separation between the coupled frequencies and offers a more accurate high-resolution PAC estimate. Moreover, the application of the pro- posed framework to the EEG data during ERN revealed significant PAC that is aligned well with prior work. In Chapter 3, we proposed a tensor based method for accurately quantifying inter-areal PAC across multiple brain regions. The proposed method employ a tensor representation of PAC across all channels, frequency bands, and subjects to quantify the dynamics of multichannel CFC. Multi- way analysis of pairwise PAC in EEG takes full advantage of the multilinear structure of the data to provide a better understanding of the interactions across multiple modes. The proposed framework of quantifying inter-areal PAC across multiple frequency bands is implemented through HoRPCA by employing the tensors constructed from PAC values across different channels, frequency bands, and subjects. The proposed method is purely data-driven and captures the spectral and spatial vari- ations in PAC, simultaneously. Also, this method offers higher accuracy and resolution compared to existing multichannel methods as it does not rely on either the Hilbert transform or the wavelet transform. The performance of the proposed measure is evaluated using simulated signals with a 108 specific application of interest in investigating whether this tensor based multichannel PAC mea- sure can successfully detect and distinguish between actual and spurious couplings. Afterward, the proposed approach was applied to multi-channel EEG recordings collected during a cognitive control study to differentiate between various response types and to identify the brain regions and frequency bands associated with those response types. Results obtained from the analyses of sim- ulated and EEG data demonstrate that our method was able to accurately filter out the background PAC activity to capture the event-related PAC dynamics. This method also detected the frequency bands responsible for these interactions while discarding the non-significant couplings. With these unique properties, the proposed multichannel PAC measure can lead to obtaining a complete un- derstanding of spatially distributed brain networks across different frequency bands. In Chapter 4, we focused on extending the proposed PAC measure to study the dynamic na- ture of PAC. Existing PAC measures are limited to quantifying the average coupling within a time window of interest. However, as PAC is dynamic, it is necessary to quantify time-varying PAC. Ex- isting time-varying PAC approaches are based on using a sliding window approach. However, this method requires the user to choose the window length and amount of overlap, which substantially affects the quantified PAC. In this chapter, In this work, we proposed a novel dynamic PAC com- putation approach that addresses these limitations. This procedure is based on matching pursuit. This approach decomposes the signal into time and frequency localized atoms that best describe the signal. Using the properties of these time and frequency localized MP atoms, we estimated the coupled low-high frequency pairs by employing k-means clustering. We quantified PAC by extracting the instantaneous phase and amplitude components from these atoms, which allowed us to capture the temporal variation of PAC using a data-driven approach. Moreover, unlike the con- ventional Hilbert transform based methods, this approach does not require any bandpass filtering. We verified the performance of the proposed method on both synthesized and real EEG data. The 109 results from synthesized data show that the proposed method detects the coupled frequencies and the time variation of the coupling correctly with high time and frequency resolution. The analysis of EEG data revealed theta-gamma and alpha-gamma PAC during response and post-response time intervals. Compared to the existing sliding window based approach, the MP based dynamic PAC measure was more effective at capturing PAC within a short time window and was more robust to noise. From these results, we posit that the proposed MP based data-driven approach offers a more effective way to quantify and track dynamic PAC. 5.2 Future Work In Chapter 2, we proposed a time frequency based PAC measure and in Chapter 3 we extended it to quantify inter-areal PAC across different brain regions and frequency bands. Although the proposed multivariate approach is promising, there are some limitations. First, the multivari- ate analysis was conducted by constructing the network with amplitude extracted only from the gamma frequency band. The low-high frequency band combinations were chosen based on the literature indicating the presence of PAC mainly between the low frequency phase (typically 5–12 Hz) and high frequency amplitude (typically 30–100Hz) [52]. Consequently, future work could include extending the proposed measure to capture other low-high frequency combinations such as delta-theta/alpha/beta, theta-alpha/beta and alpha-beta by adding those networks in A . Also the multivariate measure here focuses on the average PAC values across a pre-determined time window. However, there is growing empirical evidence that PAC is dynamic, varying across time. Consequently, future work could include extending the proposed measure by studying higher or- der tensor models to capture the temporal variation along with the spatial and spectral variation of PAC. 110 MP based dynamic PAC quantification method presented in Chapter 4 focuses on quantifying univariate PAC, i,e. PAC within a single channel. However, most neurophysiological data offers multiple channel recordings. Future research would consider incorporating the MP based dynamic PAC measure with the tensor based inter-areal PAC measure proposed in Chapter 3 to capture the dynamic PAC between spatially distributed brain regions. In this manner, we can obtain PAC across electrodes, time and, frequency. For the implementation of the MP algorithm, we have used an enriched dictionary of functions of various degrees of asymmetry – from symmetric Gabor atoms to highly asymmetric waveforms. This approach provides a flexible and sparse representation of the time series allowing us to capture any asymmetric components of the physiological signals, which are ill-described by the symmetric functions coming from the Gabor family. However, a few other dictionaries have also been proposed, for example time–frequency atoms based on damped sinusoids was shown to be more suitable for representing transients [158]. An Epsilon- Skew-Normal dictionary was proposed in [159] for the decomposition of single and multichannel EEG data. As such, in the future, the investigation of the usefulness of other dictionaries could be explored. Moreover, the multivariate matching pursuit [160, 161, 162] algorithm has been proposed to study multichannel EEG/MEG data, which could be explored in the future. Finally, the PAC considered for above experiments is non-directional. Despite the widespread use of PAC, the directionality of observed coupling remains unknown. In particular, it is not possible to determine whether the low frequency phase oscillations cause the changes in the am- plitude of the high frequency oscillations or vice versa. As neuronal interactions are directional [163], for a complete understanding of the role of PAC in neuronal transmission, estimating the directionality of the coupling is as crucial as the detection of coupling. We have observed some preliminary results based on Granger causality (GC). This approach estimates PAC directionality by quantifying the Granger causality (GC) between the phase and amplitude components extracted 111 from RID-Rihaczek t-f distribution of the coupled signals. Assuming that a driver signal oscillat- ing at a certain frequency is driving a receiver signal of another frequency and the transmission is associated with a fixed delay, the directionality of the interaction can be interpreted through GC as the driver phase/amplitude is predictive of the receiver amplitude/phase. Some preliminary experiment on simulated and real EEG data suggests that GC based approach can capture the direc- tionality of coupling from oscillatory driver to receiver during cross-frequency neural interaction. One limitation of GC based PAC directionality approach is that it can capture only linear interac- tions. However, the neural interactions in the brain is known to be nonlinear, and thus, GC based approach may fail to capture the hidden non-linear interactions in the brain. To circumvent this problem, various directed information metrics such as quantitatively measuring the PAC direction- ality by transfer entropy (TE), which has been demonstrated to be able to capture both linear and non-linear causal relationships in the brain effectively [164]. Dynamic causal modeling (DCM) could also be applied in the same manner as it has been used previously to model directed inter- actions between brain regions [165]. Future work will concentrate on the investigation of different directionality approaches to effectively infer the PAC directionality and compare their performance with GC-based PAC directionality approach. 112 BIBLIOGRAPHY 113 BIBLIOGRAPHY [1] F. L. da Silva, “Eeg and meg: relevance to neuroscience,” Neuron, vol. 80, no. 5, pp. 1112– 1128, 2013. [2] A. K. Engel, C. Gerloff, C. C. Hilgetag, and G. Nolte, “Intrinsic coupling modes: multiscale interactions in ongoing brain activity,” Neuron, vol. 80, no. 4, pp. 867–886, 2013. [3] E. L. Hall, S. E. Robson, P. G. Morris, and M. J. Brookes, “The relationship between meg and fmri,” Neuroimage, vol. 102, pp. 80–91, 2014. [4] M. L. Schölvinck, D. A. Leopold, M. J. Brookes, and P. H. Khader, “The contribution of electrophysiology to functional connectivity mapping,” Neuroimage, vol. 80, pp. 297–306, 2013. [5] P. Tewarie, M. M. Schoonheim, D. I. Schouten, C. H. Polman, L. J. Balk, B. M. Uitdehaag, J. J. Geurts, A. Hillebrand, F. Barkhof, and C. J. Stam, “Functional brain networks: linking thalamic atrophy to clinical disability in multiple sclerosis, a multimodal fmri and meg study,” Human Brain Mapping, vol. 36, no. 2, pp. 603–618, 2015. [6] C. J. Stam, “Modern network science of neurological disorders,” Nature Reviews Neuro- science, vol. 15, no. 10, p. 683, 2014. [7] E. Rodriguez, N. George, J.-P. Lachaux, J. Martinerie, B. Renault, and F. J. Varela, “Per- ception’s shadow: long-distance synchronization of human brain activity,” Nature, vol. 397, no. 6718, p. 430, 1999. [8] L. T. Trujillo, M. A. Peterson, A. W. Kaszniak, and J. J. Allen, “Eeg phase synchrony differ- ences across visual perception conditions may depend on recording and analysis methods,” Clinical Neurophysiology, vol. 116, no. 1, pp. 172–189, 2005. [9] J. M. Palva, S. Palva, and K. Kaila, “Phase synchrony among neuronal oscillations in the human cortex,” Journal of Neuroscience, vol. 25, no. 15, pp. 3962–3972, 2005. [10] P. Fries, “A mechanism for cognitive dynamics: neuronal communication through neuronal coherence,” Trends in Cognitive Sciences, vol. 9, no. 10, pp. 474–480, 2005. [11] J. E. Lisman and O. Jensen, “The theta-gamma neural code,” Neuron, vol. 77, no. 6, pp. 1002–1016, 2013. [12] R. T. Canolty, E. Edwards, S. S. Dalal, M. Soltani, S. S. Nagarajan, H. E. Kirsch, M. S. Berger, N. M. Barbaro, and R. T. Knight, “High gamma power is phase-locked to theta oscillations in human neocortex,” Science, vol. 313, no. 5793, pp. 1626–1628, 2006. 114 [13] O. Jensen and L. L. Colgin, “Cross-frequency coupling between neuronal oscillations,” Trends in Cognitive Sciences, vol. 11, no. 7, pp. 267–269, 2007. [14] R. T. Canolty and R. T. Knight, “The functional role of cross-frequency coupling,” Trends in Cognitive Sciences, vol. 14, no. 11, pp. 506–515, 2010. [15] V. Jirsa and V. Müller, “Cross-frequency coupling in real and virtual brain networks,” Fron- tiers in computational neuroscience, vol. 7, p. 78, 2013. [16] A. Hyafil, A.-L. Giraud, L. Fontolan, and B. Gutkin, “Neural cross-frequency coupling: connecting architectures, mechanisms, and functions,” Trends in Neurosciences, vol. 38, no. 11, pp. 725–740, 2015. [17] E. A. Allen, J. Liu, K. A. Kiehl, J. Gelernter, G. D. Pearlson, N. I. Perrone-Bizzozero, and V. D. Calhoun, “Components of cross-frequency modulation in health and disease,” Frontiers in systems neuroscience, vol. 5, p. 59, 2011. [18] R. van der Meij, M. Kahana, and E. Maris, “Phase–amplitude coupling in human electro- corticography is spatially distributed and phase diverse,” Journal of Neuroscience, vol. 32, no. 1, pp. 111–123, 2012. [Online]. Available: http://www.jneurosci.org/content/32/1/111 [19] M. X. Cohen, N. Axmacher, D. Lenartz, C. E. Elger, V. Sturm, and T. E. Schlaepfer, “Good vibrations: cross-frequency coupling in the human nucleus accumbens during reward pro- cessing,” Journal of Cognitive Neuroscience, vol. 21, no. 5, pp. 875–889, 2009. [20] P. Lakatos, G. Karmos, A. D. Mehta, I. Ulbert, and C. E. Schroeder, “Entrainment of neu- ronal oscillations as a mechanism of attentional selection,” Science, vol. 320, no. 5872, pp. 110–113, 2008. [21] A. B. Tort, R. W. Komorowski, J. R. Manns, N. J. Kopell, and H. Eichenbaum, “Theta– gamma coupling increases during the learning of item–context associations,” Proceedings of the National Academy of Sciences, vol. 106, no. 49, pp. 20 942–20 947, 2009. [22] J. Weule et al., “Detection of n: m phase locking from noisy data: application to magne- toencephalography,” Phys. Rev. Lett, vol. 81, no. 15, pp. 3291–3294, 1998. [23] R. T. Canolty, M. Soltani, S. S. Dalal, E. Edwards, N. F. Dronkers, S. S. Nagarajan, H. E. Kirsch, N. M. Barbaro, and R. T. Knight, “Spatiotemporal dynamics of word processing in the human brain,” Frontiers in Neuroscience, vol. 1, p. 14, 2007. [24] F. Darvas, K. J. Miller, R. P. Rao, and J. G. Ojemann, “Nonlinear phase–phase cross- frequency coupling mediates communication between distant sites in human neocortex,” Journal of Neuroscience, vol. 29, no. 2, pp. 426–435, 2009. 115 [25] A. Bruns and R. Eckhorn, “Task-related coupling from high-to low-frequency signals among visual cortical areas in human subdural recordings,” International Journal of Psychophysi- ology, vol. 51, no. 2, pp. 97–116, 2004. [26] B. Voytek, L. Secundo, A. Bidet-Caulet, D. Scabini, S. I. Stiver, A. D. Gean, G. T. Manley, and R. T. Knight, “Hemicraniectomy: a new model for human electrophysiology with high spatio-temporal resolution,” Journal of Cognitive Neuroscience, vol. 22, no. 11, pp. 2491– 2502, 2010. [27] M. Bekisz and A. Wróbel, “Coupling of beta and gamma activity in corticothalamic system of cats attending to visual stimuli,” Neuroreport, vol. 10, no. 17, pp. 3589–3594, 1999. [28] A. Bruns, R. Eckhorn, H. Jokeit, and A. Ebner, “Amplitude envelope correlation detects coupling among incoherent brain signals,” Neuroreport, vol. 11, no. 7, pp. 1509–1514, 2000. [29] M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, “Phase synchronization of chaotic oscilla- tors,” Physical review letters, vol. 76, no. 11, p. 1804, 1996. [30] F. P. De Lange, O. Jensen, M. Bauer, and I. Toni, “Interactions between posterior gamma and frontal alpha/beta oscillations during imagined actions,” Frontiers in human neuroscience, vol. 2, p. 7, 2008. [31] R. F. Helfrich, M. Huang, G. Wilson, and R. T. Knight, “Prefrontal cortex modulates pos- terior alpha oscillations during top-down guided visual perception,” Proceedings of the Na- tional Academy of Sciences, vol. 114, no. 35, pp. 9457–9462, 2017. [32] J. M. Palva and S. Palva, “Functional integration across oscillation frequencies by cross- frequency phase synchronization,” European Journal of Neuroscience, vol. 48, no. 7, pp. 2399–2406, 2018. [33] H. Witte, P. Putsche, C. Hemmelmann, C. Schelenz, and L. Leistritz, “Analysis and mod- eling of time-variant amplitude–frequency couplings of and between oscillations of eeg bursts,” Biological Cybernetics, vol. 99, no. 2, pp. 139–157, 2008. [34] S. Ray and J. H. Maunsell, “Differences in gamma frequencies across visual cortex restrict their possible use in computation,” Neuron, vol. 67, no. 5, pp. 885–896, 2010. [35] M. J. Roberts, E. Lowet, N. M. Brunet, M. Ter Wal, P. Tiesinga, P. Fries, and P. De Weerd, “Robust gamma coherence between macaque v1 and v2 by dynamic frequency matching,” Neuron, vol. 78, no. 3, pp. 523–536, 2013. [36] S. Vanhatalo, J. M. Palva, M. Holmes, J. Miller, J. Voipio, and K. Kaila, “Infras low os- cillations modulate excitability and interictal epileptic activity in the human cortex during sleep,” Proceedings of the National Academy of Sciences, vol. 101, no. 14, pp. 5053–5057, 2004. 116 [37] N. Axmacher, M. M. Henseler, O. Jensen, I. Weinreich, C. E. Elger, and J. Fell, “Cross- frequency coupling supports multi-item working memory in the human hippocampus,” Proceedings of the National Academy of Sciences, vol. 107, no. 7, pp. 3228–3233, 2010. [Online]. Available: http://www.pnas.org/content/107/7/3228 [38] G. Buzsáki, D. Buhl, K. Harris, J. Csicsvari, B. Czeh, and A. Morozov, “Hippocampal network patterns of activity in the mouse,” Neuroscience, vol. 116, no. 1, pp. 201–211, 2003. [39] H. Hentschke, M. G. Perkins, R. A. Pearce, and M. I. Banks, “Muscarinic blockade weakens interaction of gamma with theta rhythms in mouse hippocampus,” European Journal of Neuroscience, vol. 26, no. 6, pp. 1642–1656, 2007. [40] A. Bragin, G. Jandó, Z. Nádasdy, J. Hetke, K. Wise, and G. Buzsáki, “Gamma (40-100 hz) oscillation in the hippocampus of the behaving rat,” Journal of Neuroscience, vol. 15, no. 1, pp. 47–60, 1995. [41] K. M. Kendrick, Y. Zhan, H. Fischer, A. U. Nicol, X. Zhang, and J. Feng, “Learning alters theta amplitude, theta-gamma coupling and neuronal synchronization in inferotemporal cortex,” BMC Neuroscience, vol. 12, no. 1, p. 55, Jun 2011. [Online]. Available: https://doi.org/10.1186/1471-2202-12-55 [42] P. Lakatos, A. S. Shah, K. H. Knuth, I. Ulbert, G. Karmos, and C. E. Schroeder, “An oscillatory hierarchy controlling neuronal excitability and stimulus processing in the auditory cortex,” Journal of Neurophysiology, vol. 94, no. 3, pp. 1904–1911, 2005, pMID: 15901760. [Online]. Available: https://doi.org/10.1152/jn.00263.2005 [43] M. X. Cohen, C. E. Elger, and J. Fell, “Oscillatory activity and phase–amplitude coupling in the human medial frontal cortex during decision making,” Journal of Cognitive Neuro- science, vol. 21, no. 2, pp. 390–402, 2008. [44] C. E. Schroeder and P. Lakatos, “Low-frequency neuronal oscillations as instruments of sensory selection,” Trends in Neurosciences, vol. 32, no. 1, pp. 9 – 18, 2009. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0166223608002506 [45] B. Händel and T. Haarmeier, “Cross-frequency coupling of brain oscillations indicates the success in visual motion discrimination,” Neuroimage, vol. 45, no. 3, pp. 1040–1046, 2009. [46] T. Demiralp, Z. Bayraktaroglu, D. Lenz, S. Junge, N. A. Busch, B. Maess, M. Ergen, and C. S. Herrmann, “Gamma amplitudes are coupled to theta phase in human eeg during visual perception,” International Journal of Psychophysiology, vol. 64, no. 1, pp. 24 – 30, 2007, brain Oscillations:Cutting Edges. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S016787600600211X 117 [47] S. Khan, A. Gramfort, N. R. Shetty, M. G. Kitzbichler, S. Ganesan, J. M. Moran, S. M. Lee, J. D. E. Gabrieli, H. B. Tager-Flusberg, R. M. Joseph, M. R. Herbert, M. S. Hämäläinen, and T. Kenet, “Local and long-range functional connectivity is reduced in concert in autism spectrum disorders,” Proceedings of the National Academy of Sciences, vol. 110, no. 8, pp. 3107–3112, 2013. [Online]. Available: http://www.pnas.org/content/110/8/3107 [48] K. Kirihara, A. J. Rissling, N. R. Swerdlow, D. L. Braff, and G. A. Light, “Hierarchical organization of gamma and theta oscillatory dynamics in schizophrenia,” Biological Psychi- atry, vol. 71, no. 10, pp. 873 – 880, 2012, altered Functional Connectivity in Schizophrenia. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0006322312000558 [49] C. de Hemptinne, E. S. Ryapolova-Webb, E. L. Air, P. A. Garcia, K. J. Miller, J. G. Ojemann, J. L. Ostrem, N. B. Galifianakis, and P. A. Starr, “Exaggerated phase–amplitude coupling in the primary motor cortex in parkinson disease,” Proceedings of the National Academy of Sciences, vol. 110, no. 12, pp. 4780–4785, 2013. [Online]. Available: http://www.pnas.org/content/110/12/4780 [50] R. F. Helfrich and R. T. Knight, “Oscillatory dynamics of prefrontal cognitive control,” Trends in cognitive sciences, vol. 20, no. 12, pp. 916–930, 2016. [51] J. Riddle, D. A. Vogelsang, K. Hwang, D. Cellier, and M. D’Esposito, “Distinct oscilla- tory dynamics underlie different components of hierarchical cognitive control,” Journal of Neuroscience, vol. 40, no. 25, pp. 4945–4953, 2020. [52] A. B. Tort, R. Komorowski, H. Eichenbaum, and N. Kopell, “Measuring phase-amplitude coupling between neuronal oscillations of different frequencies,” Journal of Neurophysiol- ogy, vol. 104, no. 2, pp. 1195–1210, 2010. [53] T. Dupré la Tour, L. Tallot, L. Grabot, V. Doyère, V. van Wassenhove, Y. Grenier, and A. Gramfort, “Non-linear auto-regressive models for cross-frequency coupling in neural time series,” PLOS Computational Biology, vol. 13, no. 12, pp. 1–32, 12 2017. [Online]. Available: https://doi.org/10.1371/journal.pcbi.1005893 [54] M. J. Huelsemann, E. Naumann, and B. Rasch, “Quantification of phase- amplitude coupling in neuronal oscillations: Comparison of phase-locking value, mean vector length, and modulation index,” BioRxiv, 2018. [Online]. Available: https://www.biorxiv.org/content/early/2018/03/28/290361 [55] J. Aru, J. Aru, V. Priesemann, M. Wibral, L. Lana, G. Pipa, W. Singer, and R. Vicente, “Untangling cross-frequency coupling in neuroscience,” Current Opinion in Neurobiology, vol. 31, pp. 51 – 61, 2015, sI: Brain rhythms and dynamic coordination. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0959438814001640 118 [56] B. van Wijk, A. Jha, W. Penny, and V. Litvak, “Parametric estimation of cross-frequency coupling,” Journal of Neuroscience Methods, vol. 243, pp. 94 – 102, 2015. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0165027015000394 [57] M. X. Cohen, “Multivariate cross-frequency coupling via generalized eigendecomposition,” eLife, vol. 6, p. e21792, jan 2017. [Online]. Available: https://doi.org/10.7554/eLife.21792 [58] D. Dvorak and A. A. Fenton, “Toward a proper estimation of phase–amplitude coupling in neural oscillations,” Journal of Neuroscience Methods, vol. 225, pp. 42 – 56, 2014. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0165027014000132 [59] A. Hyafil, “Misidentifications of specific forms of cross-frequency coupling: three warnings,” Frontiers in Neuroscience, vol. 9, p. 370, 2015. [Online]. Available: https://www.frontiersin.org/article/10.3389/fnins.2015.00370 [60] M. A. Kramer, A. B. Tort, and N. J. Kopell, “Sharp edge artifacts and spurious coupling in eeg frequency comodulation measures,” Journal of Neuro- science Methods, vol. 170, no. 2, pp. 352 – 357, 2008. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0165027008000538 [61] J. I. Berman, J. McDaniel, S. Liu, L. Cornew, W. Gaetz, T. P. Roberts, and J. C. Edgar, “Variable bandwidth filtering for improved sensitivity of cross-frequency coupling metrics,” Brain Connectivity, vol. 2, no. 3, pp. 155–163, 2012, pMID: 22577870. [Online]. Available: https://doi.org/10.1089/brain.2012.0085 [62] A. Nakhnikian, S. Ito, L. Dwiel, L. Grasse, G. Rebec, L. Lauridsen, and J. Beggs, “A novel cross-frequency coupling detection method using the generalized morse wavelets,” Journal of neuroscience methods, vol. 269, pp. 61–73, 2016. [63] Y. Zerouali, J.-M. Lina, Z. Sekerovic, J. Godbout, J. Dube, P. Jolicoeur, and J. Carrier, “A time-frequency analysis of the dynamics of cortical networks of sleep spindles from meg- eeg recordings,” Frontiers in neuroscience, vol. 8, p. 310, 2014. [64] T. E. Özkurt and A. Schnitzler, “A critical note on the definition of phase–amplitude cross- frequency coupling,” Journal of Neuroscience methods, vol. 201, no. 2, pp. 438–443, 2011. [65] M. X. Cohen, “Assessing transient cross-frequency coupling in eeg data,” Journal of neuro- science methods, vol. 168, no. 2, pp. 494–499, 2008. [66] F. Mormann, J. Fell, N. Axmacher, B. Weber, K. Lehnertz, C. E. Elger, and G. Fernández, “Phase/amplitude reset and theta–gamma interaction in the human medial temporal lobe during a continuous word recognition memory task,” Hippocampus, vol. 15, no. 7, pp. 890– 900, 2011. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/hipo.20117 119 [67] A. Rihaczek, “Signal energy distribution in time and frequency,” IEEE Transactions on Information Theory, vol. 14, no. 3, pp. 369–374, May 1968. [68] S. Aviyente and A. Y. Mutlu, “A time-frequency-based approach to phase and phase syn- chrony estimation,” IEEE Transactions on Signal Processing, vol. 59, no. 7, pp. 3086–3098, July 2011. [69] S. Aviyente, E. M. Bernat, W. S. Evans, and S. R. Sponheim, “A phase synchrony measure for quantifying dynamic functional integration in the brain,” Human Brain Mapping, vol. 32, no. 1, pp. 80–93, 2011. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/hbm.21000 [70] G. Buzsaki, Rhythms of the Brain. Oxford University Press, 2006. [71] C. De Hemptinne, E. S. Ryapolova-Webb, E. L. Air, P. A. Garcia, K. J. Miller, J. G. Oje- mann, J. L. Ostrem, N. B. Galifianakis, and P. A. Starr, “Exaggerated phase–amplitude coupling in the primary motor cortex in parkinson disease,” Proceedings of the National Academy of Sciences, p. 201214546, 2013. [72] G. Buzsáki, “Neural syntax: cell assemblies, synapsembles, and readers,” Neuron, vol. 68, no. 3, pp. 362–385, 2010. [73] M. A. Kramer, A. B. Tort, and N. J. Kopell, “Sharp edge artifacts and spurious coupling in eeg frequency comodulation measures,” Journal of neuroscience methods, vol. 170, no. 2, pp. 352–357, 2008. [74] A. Hyafil, “Misidentifications of specific forms of cross-frequency coupling: three warn- ings,” Frontiers in neuroscience, vol. 9, p. 370, 2015. [75] M. X. Cohen, “Multivariate cross-frequency coupling via generalized eigendecomposition,” ELife, vol. 6, p. e21792, 2017. [76] S. Aviyente and A. Y. Mutlu, “A time-frequency-based approach to phase and phase syn- chrony estimation,” IEEE Transactions on Signal Processing, vol. 59, no. 7, pp. 3086–3098, 2011. [77] B. Voytek, M. D’Esposito, N. Crone, and R. T. Knight, “A method for event-related phase/amplitude coupling,” NeuroImage, vol. 64, pp. 416 – 424, 2013. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S1053811912009330 [78] B. Pittman-Polletta, W.-H. Hsieh, S. Kaur, M.-T. Lo, and K. Hu, “Detecting phase-amplitude coupling with high frequency resolution using adaptive decompositions,” Journal of neuro- science methods, vol. 226, pp. 15–32, 2014. 120 [79] T. P. Moran, E. M. Bernat, S. Aviyente, H. S. Schroder, and J. S. Moser, “Sending mixed signals: worry is associated with enhanced initial error processing but reduced call for sub- sequent cognitive control,” Social Cognitive and Affective neuroscience, vol. 10, no. 11, pp. 1548–1556, 2015. [80] S. Samiee and S. Baillet, “Time-resolved phase-amplitude coupling in neural oscillations,” NeuroImage, vol. 159, pp. 270–279, 2017. [81] K. J. Miller, L. B. Sorensen, J. G. Ojemann, and M. den Nijs, “Power-law scaling in the brain surface electric potential,” PLOS Computational Biology, vol. 5, no. 12, pp. 1–10, 12 2009. [Online]. Available: https://doi.org/10.1371/journal.pcbi.1000609 [82] B. J. He, J. M. Zempel, A. Z. Snyder, and M. E. Raichle, “The temporal structures and func- tional significance of scale-free brain activity,” Neuron, vol. 66, no. 3, pp. 353 – 369, 2010. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0896627310002916 [83] H. Zhivomirov, “A method for colored noise generation,” Romanian Journal of Acoustics and Vibration, vol. 15, no. 1, pp. 14–19, 2018. [84] M. Kramer and U. Eden, “Assessment of cross-frequency coupling with confidence using generalized linear models,” Journal of Neuroscience Meth- ods, vol. 220, no. 1, pp. 64 – 74, 2013. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S016502701300277X [85] M. L. V. Quyen, J. Foucher, J.-P. Lachaux, E. Rodriguez, A. Lutz, J. Martinerie, and F. J. Varela, “Comparison of hilbert transform and wavelet methods for the analysis of neuronal synchrony,” Journal of Neuroscience Methods, vol. 111, no. 2, pp. 83 – 98, 2001. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0165027001003727 [86] R. A. Seymour, G. Rippon, and K. Kessler, “The detection of phase amplitude coupling during sensory processing,” Frontiers in Neuroscience, vol. 11, p. 487, 2017. [Online]. Available: https://www.frontiersin.org/article/10.3389/fnins.2017.00487 [87] L. Stanković, “A measure of some time–frequency distributions concentration,” Signal Pro- cessing, vol. 81, no. 3, pp. 621–631, 2001. [88] R. G. Baraniuk, P. Flandrin, A. J. Janssen, and O. J. Michel, “Measuring time-frequency information content using the rényi entropies,” IEEE Transactions on Information Theory, vol. 47, no. 4, pp. 1391–1409, 2001. [89] B. A. Eriksen and C. W. Eriksen, “Effects of noise letters upon the identification of a target letter in a nonsearch task,” Perception & Psychophysics, vol. 16, no. 1, pp. 143–149, Jan 1974. [Online]. Available: https://doi.org/10.3758/BF03203267 121 [90] J. Kayser and C. E. Tenke, “Principal components analysis of laplacian waveforms as a generic method for identifying erp generator patterns: I. evaluation with auditory oddball tasks,” Clinical neurophysiology, vol. 117, no. 2, pp. 348–368, 2006. [91] B. Voytek, R. Canolty, A. Shestyuk, N. Crone, J. Parvizi, and R. Knight, “Shifts in gamma phase–amplitude coupling frequency from theta to alpha over posterior cortex during visual tasks,” Frontiers in Human Neuroscience, vol. 4, p. 191, 2010. [Online]. Available: https://www.frontiersin.org/article/10.3389/fnhum.2010.00191 [92] M. Bonnefond and O. Jensen, “Gamma activity coupled to alpha phase as a mechanism for top-down controlled gating,” PLOS ONE, vol. 10, no. 6, pp. 1–11, 06 2015. [Online]. Available: https://doi.org/10.1371/journal.pone.0128667 [93] E. Spaak, M. Bonnefond, A. Maier, D. A. Leopold, and O. Jensen, “Layer-specific entrainment of gamma-band neural activity by the alpha rhythm in monkey visual cortex,” Current Biology, vol. 22, no. 24, pp. 2313 – 2318, 2012. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0960982212012055 [94] M. X. Cohen and S. van Gaal, “Dynamic interactions between large-scale brain networks predict behavioral adaptation after perceptual errors,” Cerebral Cortex, vol. 23, no. 5, pp. 1061–1072, 2013. [Online]. Available: http://dx.doi.org/10.1093/cercor/bhs069 [95] L. T. Trujillo and J. J. Allen, “Theta eeg dynamics of the error-related negativity,” Clinical Neurophysiology, vol. 118, no. 3, pp. 645 – 668, 2007. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S1388245706015264 [96] C. B. Holroyd and M. G. Coles, “The neural basis of human error processing: reinforcement learning, dopamine, and the error-related negativity.” Psychological Review, vol. 109, no. 4, p. 679, 2002. [97] K. R. Ridderinkhof, M. Ullsperger, E. A. Crone, and S. Nieuwenhuis, “The role of the medial frontal cortex in cognitive control,” Science, vol. 306, no. 5695, pp. 443–447, 2004. [98] M. M. Botvinick, T. S. Braver, D. M. Barch, C. S. Carter, and J. D. Cohen, “Conflict moni- toring and cognitive control.” Psychological Review, vol. 108, no. 3, p. 624, 2001. [99] S. Aviyente and W. J. Williams, “A centrosymmetric kernel decomposition for time- frequency distribution computation,” IEEE Transactions on Signal Processing, vol. 52, no. 6, pp. 1574–1584, 2004. [100] A. de Cheveigné and L. C. Parra, “Joint decorrelation, a versatile tool for multichannel data analysis,” NeuroImage, vol. 98, pp. 487 – 505, 2014. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S1053811914004534 122 [101] P. Tewarie, A. Hillebrand, B. W. van Dijk, C. J. Stam, G. C. O’Neill, P. V. Mieghem, J. M. Meier, M. W. Woolrich, P. G. Morris, and M. J. Brookes, “Integrating cross-frequency and within band functional networks in resting-state meg: A multi-layer network approach,” NeuroImage, vol. 142, pp. 324 – 336, 2016. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S1053811916303718 [102] M. J. Brookes, P. K. Tewarie, B. A. Hunt, S. E. Robson, L. E. Gascoyne, E. B. Liddle, P. F. Liddle, and P. G. Morris, “A multi-layer network approach to meg connectivity analysis,” NeuroImage, vol. 132, pp. 425 – 438, 2016. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S1053811916001543 [103] A. B. Tort, M. A. Kramer, C. Thorn, D. J. Gibson, Y. Kubota, A. M. Graybiel, and N. J. Kopell, “Dynamic cross-frequency couplings of local field potential oscillations in rat stria- tum and hippocampus during performance of a t-maze task,” Proceedings of the National Academy of Sciences, vol. 105, no. 51, pp. 20 517–20 522, 2008. [104] N. Axmacher, M. M. Henseler, O. Jensen, I. Weinreich, C. E. Elger, and J. Fell, “Cross- frequency coupling supports multi-item working memory in the human hippocampus,” Pro- ceedings of the National Academy of Sciences, vol. 107, no. 7, pp. 3228–3233, 2010. [105] R. T. Canolty, C. F. Cadieu, K. Koepsell, R. T. Knight, and J. M. Carmena, “Multivariate phase–amplitude cross-frequency coupling in neurophysiological signals,” IEEE Transac- tions on Biomedical Engineering, vol. 59, no. 1, pp. 8–11, 2011. [106] T. T. K. Munia and S. Aviyente, “A time-frequency based multivariate phase-amplitude coupling measure,” in ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), May 2019, pp. 1095–1099. [107] Y.-P. Lin, P.-K. Jao, and Y.-H. Yang, “Improving cross-day eeg-based emotion classifica- tion using robust principal component analysis,” Frontiers in computational neuroscience, vol. 11, p. 64, 2017. [108] C.-S. Wei, Y.-P. Lin, and T.-P. Jung, “Exploring the eeg correlates of neurocognitive lapse with robust principal component analysis,” in International Conference on Augmented Cog- nition. Springer, 2016, pp. 113–120. [109] A. Ozdemir, E. M. Bernat, and S. Aviyente, “Recursive tensor subspace tracking for dy- namic brain network analysis,” IEEE Transactions on Signal and Information Processing over Networks, vol. 3, no. 4, pp. 669–682, 2017. [110] J. Daume, T. Gruber, A. K. Engel, and U. Friese, “Phase-amplitude coupling and long-range phase synchronization reveal frontotemporal interactions during visual working memory,” Journal of Neuroscience, vol. 37, no. 2, pp. 313–322, 2017. 123 [111] M. X. Cohen and S. Van Gaal, “Dynamic interactions between large-scale brain networks predict behavioral adaptation after perceptual errors,” Cerebral Cortex, vol. 23, no. 5, pp. 1061–1072, 2013. [112] A. G. Mahyari and S. Aviyente, “Hierarchical distributed compressive sensing for simulta- neous sparse approximation and common component extraction,” Digital Signal Processing, vol. 60, pp. 230–241, 2017. [113] “https://github.com/muntam/tf-pac.” 2019. [Online]. Available: https://github.com/muntam/TF-PAC [114] J. Aru, J. Aru, V. Priesemann, M. Wibral, L. Lana, G. Pipa, W. Singer, and R. Vicente, “Untangling cross-frequency coupling in neuroscience,” Current opinion in neurobiology, vol. 31, pp. 51–61, 2015. [115] F. Miwakeichi, E. Martınez-Montes, P. A. Valdés-Sosa, N. Nishiyama, H. Mizuhara, and Y. Yamaguchi, “Decomposing eeg data into space–time–frequency components using par- allel factor analysis,” NeuroImage, vol. 22, no. 3, pp. 1035–1045, 2004. [116] H. Lee, Y.-D. Kim, A. Cichocki, and S. Choi, “Nonnegative tensor factorization for continu- ous eeg classification,” International journal of neural systems, vol. 17, no. 04, pp. 305–317, 2007. [117] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” SIAM review, vol. 51, no. 3, pp. 455–500, 2009. [118] M. E. Timmerman and H. A. Kiers, “Three-mode principal components analysis: Choosing the numbers of components and sensitivity to local optima,” British journal of mathematical and statistical psychology, vol. 53, no. 1, pp. 1–16, 2000. [119] F. Cong, Q.-H. Lin, L.-D. Kuang, X.-F. Gong, P. Astikainen, and T. Ristaniemi, “Tensor decomposition of eeg signals: a brief review,” Journal of neuroscience methods, vol. 248, pp. 59–69, 2015. [120] E. J. Candès, X. Li, Y. Ma, and J. Wright, “Robust principal component analysis?” Journal of the ACM (JACM), vol. 58, no. 3, pp. 1–37, 2011. [121] D. Goldfarb and Z. Qin, “Robust low-rank tensor recovery: Models and algorithms,” SIAM Journal on Matrix Analysis and Applications, vol. 35, no. 1, pp. 225–253, 2014. [122] J. Li, G. Han, J. Wen, and X. Gao, “Robust tensor subspace learning for anomaly detection,” International Journal of Machine Learning and Cybernetics, vol. 2, no. 2, pp. 89–98, 2011. [123] H. Tan, J. Feng, G. Feng, W. Wang, and Y.-J. Zhang, “Traffic volume data outlier recovery via tensor model,” Mathematical Problems in Engineering, vol. 2013, 2013. 124 [124] K. Xie, X. Li, X. Wang, G. Xie, J. Wen, J. Cao, and D. Zhang, “Fast tensor factorization for accurate internet anomaly detection,” IEEE/ACM transactions on networking, vol. 25, no. 6, pp. 3794–3807, 2017. [125] M. J. Hülsemann, E. Naumann, and B. Rasch, “Quantification of phase-amplitude coupling in neuronal oscillations: Comparison of phase-locking value, mean vector length, modula- tion index, and generalized linear modeling cross-frequency coupling,” Frontiers in neuro- science, vol. 13, p. 573, 2019. [126] A. Gramfort, T. Papadopoulo, E. Olivi, and M. Clerc, “Openmeeg: opensource software for quasistatic bioelectromagnetics,” Biomedical engineering online, vol. 9, no. 1, p. 45, 2010. [127] F. Tadel, S. Baillet, J. C. Mosher, D. Pantazis, and R. M. Leahy, “Brainstorm: a user-friendly application for meg/eeg analysis,” Computational intelligence and neuroscience, vol. 2011, p. 8, 2011. [128] B. A. Eriksen and C. W. Eriksen, “Effects of noise letters upon the identification of a target letter in a nonsearch task,” Perception & psychophysics, vol. 16, no. 1, pp. 143–149, 1974. [129] P. L. Nunez, R. Srinivasan, A. F. Westdorp, R. S. Wijesinghe, D. M. Tucker, R. B. Silber- stein, and P. J. Cadusch, “Eeg coherency: I: statistics, reference electrode, volume conduc- tion, laplacians, cortical imaging, and interpretation at multiple scales,” Electroencephalog- raphy and clinical neurophysiology, vol. 103, no. 5, pp. 499–515, 1997. [130] J. Wright, A. Y. Yang, A. Ganesh, S. S. Sastry, and Y. Ma, “Robust face recognition via sparse representation,” IEEE transactions on pattern analysis and machine intelligence, vol. 31, no. 2, pp. 210–227, 2008. [131] X. Luan, B. Fang, L. Liu, W. Yang, and J. Qian, “Extracting sparse error of robust pca for face recognition in the presence of varying illumination and occlusion,” Pattern Recognition, vol. 47, no. 2, pp. 495–508, 2014. [132] A. K. Engel, P. Fries, and W. Singer, “Dynamic predictions: oscillations and synchrony in top–down processing,” Nature Reviews Neuroscience, vol. 2, no. 10, pp. 704–716, 2001. [133] G. Buzsáki and A. Draguhn, “Neuronal oscillations in cortical networks,” science, vol. 304, no. 5679, pp. 1926–1929, 2004. [134] R. G. Port, M. A. Dipiero, M. Ku, S. Liu, L. Blaskey, E. S. Kuschner, J. C. Edgar, T. P. Roberts, and J. I. Berman, “Children with autism spectrum disorder demonstrate regionally specific altered resting-state phase–amplitude coupling,” Brain connectivity, vol. 9, no. 5, pp. 425–436, 2019. [135] K.-m. An, T. Ikeda, C. Hasegawa, Y. Yoshimura, S. Tanaka, D. N. Saito, K. Yaoi, S. Iwasaki, T. Hirosawa, O. Jensen et al., “Aberrant brain oscillatory coupling from the primary mo- 125 tor cortex in children with autism spectrum disorders,” NeuroImage: Clinical, vol. 29, p. 102560, 2021. [136] T. H. Lee, M. Kim, W. J. Hwang, T. Kim, Y. B. Kwak, and J. S. Kwon, “Relationship be- tween resting-state theta phase-gamma amplitude coupling and neurocognitive functioning in patients with first-episode psychosis,” Schizophrenia research, vol. 216, pp. 154–160, 2020. [137] R. Gong, M. Wegscheider, C. Mühlberg, R. Gast, C. Fricke, J.-J. Rumpf, V. V. Nikulin, T. R. Knösche, and J. Classen, “Spatiotemporal features of β -γ phase-amplitude coupling in parkinson’s disease derived from scalp eeg,” Brain, vol. 144, no. 2, pp. 487–503, 2021. [138] F. Mormann, J. Fell, N. Axmacher, B. Weber, K. Lehnertz, C. E. Elger, and G. Fernández, “Phase/amplitude reset and theta–gamma interaction in the human medial temporal lobe during a continuous word recognition memory task,” Hippocampus, vol. 15, no. 7, pp. 890– 900, 2005. [139] T. D. La Tour, L. Tallot, L. Grabot, V. Doyère, V. Van Wassenhove, Y. Grenier, and A. Gram- fort, “Non-linear auto-regressive models for cross-frequency coupling in neural time series,” PLoS computational biology, vol. 13, no. 12, p. e1005893, 2017. [140] C. S. Zandvoort and G. Nolte, “Defining the filter parameters for phase-amplitude coupling from a bispectral point of view,” Journal of Neuroscience Methods, vol. 350, p. 109032, 2021. [141] K. K. Sreenivasan, C. E. Curtis, and M. D’Esposito, “Revisiting the role of persistent neural activity during working memory,” Trends in cognitive sciences, vol. 18, no. 2, pp. 82–89, 2014. [142] K. J. Friston, “The labile brain. i. neuronal transients and nonlinear coupling,” Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences, vol. 355, no. 1394, pp. 215–236, 2000. [143] W. Penny, E. Duzel, K. Miller, and J. Ojemann, “Testing for nested oscillation,” Journal of neuroscience methods, vol. 174, no. 1, pp. 50–61, 2008. [144] B. Voytek, M. D’esposito, N. Crone, and R. T. Knight, “A method for event-related phase/amplitude coupling,” Neuroimage, vol. 64, pp. 416–424, 2013. [145] T. Spustek, W. W. Jedrzejczak, and K. J. Blinowska, “Matching pursuit with asymmet- ric functions for signal decomposition and parameterization,” PloS one, vol. 10, no. 6, p. e0131007, 2015. [146] E. Bedrosian, “A product theorem for hilbert transforms,” Proceedings of the IEEE, vol. 51, no. 5, pp. 868–869, 1963. 126 [147] S. G. Mallat and Z. Zhang, “Matching pursuits with time-frequency dictionaries,” IEEE Transactions on signal processing, vol. 41, no. 12, pp. 3397–3415, 1993. [148] A. Delorme and S. Makeig, “Eeglab: an open source toolbox for analysis of single-trial eeg dynamics including independent component analysis,” Journal of neuroscience methods, vol. 134, no. 1, pp. 9–21, 2004. [149] S. Aviyente, E. M. Bernat, S. M. Malone, and W. G. Iacono, “Time-frequency data reduction for event related potentials: combining principal component analysis and matching pursuit,” EURASIP journal on advances in signal processing, vol. 2010, p. 4, 2010. [150] S. Aviyente, “Compressed sensing framework for eeg compression,” in Proc. IEEE/SP 14th Workshop Stat. Signal Process, 2007, pp. 181–184. [151] S. Lloyd, “Least squares quantization in pcm,” IEEE transactions on information theory, vol. 28, no. 2, pp. 129–137, 1982. [152] K. D. Joshi and P. Nalwade, “Modified k-means for better initial cluster centres,” Inter- national Journal of Computer Science and Mobile Computing, vol. 2, no. 7, pp. 219–223, 2013. [153] J. I. Berman, J. McDaniel, S. Liu, L. Cornew, W. Gaetz, T. P. Roberts, and J. C. Edgar, “Variable bandwidth filtering for improved sensitivity of cross-frequency coupling metrics,” Brain connectivity, vol. 2, no. 3, pp. 155–163, 2012. [154] B. van Wijk, A. Jha, W. Penny, and V. Litvak, “Parametric estimation of cross-frequency coupling,” Journal of neuroscience methods, vol. 243, pp. 94–102, 2015. [155] S. Aviyente, E. M. Bernat, W. S. Evans, and S. R. Sponheim, “A phase synchrony measure for quantifying dynamic functional integration in the brain,” Wiley Online Library, Tech. Rep. 1, 2011. [156] T. T. Munia and S. Aviyente, “Time-frequency based phase-amplitude coupling measure for neuronal oscillations,” Scientific reports, vol. 9, no. 1, pp. 1–15, 2019. [157] T. T. K. Munia and S. Aviyente, “Multivariate analysis of bivariate phase-amplitude cou- pling in eeg data using tensor robust pca,” IEEE Transactions on Neural Systems and Reha- bilitation Engineering, vol. 29, pp. 1268–1279, 2021. [158] S. Ghofrani, D. C. McLernon, and A. Ayatollahi, “On conditional spectral moments of gaus- sian and damped sinusoidal atoms in adaptive signal decomposition,” Signal processing, vol. 85, no. 10, pp. 1984–1992, 2005. [159] D. Strohmeier, A. Halbleib, M. Gratkowski, and J. Haueisen, “The epsilon-skew-normal dictionary for the decomposition of single-and multichannel biomedical recordings using 127 matching pursuit algorithms,” in XII Mediterranean Conference on Medical and Biological Engineering and Computing 2010. Springer, 2010, pp. 97–100. [160] C. Sieluzycki, R. Konig, A. Matysiak, R. Kus, D. Ircha, and P. J. Durka, “Single-trial evoked brain responses modeled by multivariate matching pursuit,” IEEE Transactions on Biomed- ical Engineering, vol. 56, no. 1, pp. 74–82, 2008. [161] R. Kuś, P. T. Różański, and P. J. Durka, “Multivariate matching pursuit in optimal gabor dictionaries: theory and software with interface for eeg/meg via svarog,” Biomedical engi- neering online, vol. 12, no. 1, pp. 1–28, 2013. [162] R. Liu, B. Karumuri, J. Adkinson, T. N. Hutson, I. Vlachos, and L. Iasemidis, “Multivariate matching pursuit decomposition and normalized gabor entropy for quantification of preictal trends in epilepsy,” Entropy, vol. 20, no. 6, p. 419, 2018. [163] B. Nandi, P. Swiatek, B. Kocsis, and M. Ding, “Inferring the direction of rhythmic neural transmission via inter-regional phase-amplitude coupling (ir-pac),” Scientific reports, vol. 9, no. 1, pp. 1–13, 2019. [164] R. Vicente, M. Wibral, M. Lindner, and G. Pipa, “Transfer entropy—a model-free measure of effective connectivity for the neurosciences,” Journal of computational neuroscience, vol. 30, no. 1, pp. 45–67, 2011. [165] B. C. van Wijk, V. Litvak, K. J. Friston, and A. Daffertshofer, “Nonlinear coupling between occipital and motor cortex during motor imagery: a dynamic causal modeling study,” Neu- roimage, vol. 71, pp. 104–113, 2013. 128