Cortico-cortical and Cortico-muscular Connectivity Analysis: Methods and Application to Parkinson’s Disease by Joyce Hsien-yin Chiang B.A.Sc., The University of British Columbia, 2004 M.A.Sc., The University of British Columbia, 2007 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in THE FACULTY OF GRADUATE STUDIES (Electrical and Computer Engineering) The University of British Columbia (Vancouver) April 2012 c© Joyce Hsien-yin Chiang, 2012 Abstract The concept of brain connectivity provides a new perspective to the understanding of the mechanism underlying brain functions, complementing the traditional approach of analyzing neural activity of isolated regions. Among the existing connectivity analysis techniques, mul- tivariate autoregressive (mAR)-based measures are of great interest for their ability to char- acterize both directionality and spectral property of cortical interactions. Yet, the direct es- timation of mAR-based connectivity from scalp electroencephalogram (EEG) is confounded by volume conduction, statistical instability and inter-subject variability. In this thesis, we propose novel signal processing methods to enhance the existing mAR-based connectivity methods. First, we explore incorporating sparsity constraints into the mAR formulation at both subject level and group level using LASSO-based regression. We show by simulation that sparse mAR yields more stable and accurate connectivity estimates compared to the tra- ditional, non-sparse approach. Furthermore, the group-wise sparsity simplifies the inference of group-level connectivity patterns from multi-subject data. To mitigate the effect of volume conduction, we investigate source-level connectivity and propose a state-space generalized mAR framework to jointly model the mixing effect of volume conduction and causal rela- tionships between underlying neural sources. By jointly estimating the mixing process and mAR model parameters, the proposed technique demonstrates improved connectivity estima- tion performance. Finally, we expanded our connectivity analysis to cortico-muscular level by modeling the relationships between EEG and simultaneously recorded electromyography (EMG) data using a multiblock partial least square (mbPLS) framework. The hierarchical construction of the mbPLS framework provides a natural way to model multi-subject, multi- ii modal data, enabling the identification of maximally covarying common patterns from EEG and EMG across subjects. Applications of the proposed techniques to EEG and EMG data of healthy and Parkinson’s disease (PD) subjects demonstrate that directional connectivity analysis is a more sensitive technique than traditional univariate spectral analysis in revealing complex effects of motor tasks and disease. Moreover, alternations in connectivity accurately predict disease severity in PD. These new analysis tools allow a better understanding of brain function and provide a basis for developing objective measures to assess progression of neu- rological diseases. iii Preface This thesis documents the research work performed under the supervision of Dr. Z. Jane Wang. Each of Chapters 2 to 7 is based on manuscripts that have been accepted, submitted, or to be submitted for publication in international peer-reviewed journals and conferences. A version of Chapter 2 has been published as: J. Chiang, Z. J. Wang, and M. J. McKeown, “Sparse multivariate autoregressive (mAR)-based partial directed coherence (PDC) for elec- troencephalogram (EEG) analysis,” in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP’09), Apr. 2009, pp. 457-460. I was responsible for the development of analysis techniques, numerical simulation and analysis and compilation of results. I was involved in the collection of EEG data together with Magdalena Dybalski. I prepared the manuscript with subsequent editorial input from Dr. Z. Jane Wang and Dr. Martin J. McKeown. A version of Chapter 3 has been published as: G. Tropini, J. Chiang, Z. J. Wang, E. Ty, and M. J. McKeown, “Altered directional connectivity in Parkinson’s disease during performance of a visually guided task,” NeuroImage, vol. 56, no. 4, pp. 2144-2156, Jun. 2011. I was responsible for developing analysis tools and implementing trigger pulses for synchronizing visual stimulus generation with data collection. The analysis of EEG data and preparation of the manuscript were performed jointly with Giorgia Tropini with the guidance of Dr. Martin J. McKeown. A version of Chapter 4 has been submitted for publication. I was responsible for the development of analysis techniques, analysis of EEG data and preparation of the manuscript. The implementation of the experimental paradigm and collection of EEG data in Chapters 3 iv and 4 were performed by Giorgia Tropini. A version of Chapter 5 has been submitted for publication. The design of the experiment and data collection were carried out by Vignan Yogendrakumar. I was responsible for the data preprocessing, development of analysis techniques and analysis of results. The manuscript was drafted jointly with Vignan Yogendrakumar and Diana Kim with the guidance of Dr. Martin J. McKeown. A version of Chapter 6 has been published as: J. Chiang, Z. J. Wang, and M. J. McKeown, “A generalized multivariate autoregressive (GmAR)-based approach for EEG source connec- tivity analysis,” IEEE Transactions on Signal Processing, vol. 60, no. 1, pp. 453–465, Jan. 2012. Part of the content in Chapter 7 has been published as: J. Chiang, Z. J. Wang, and M. J. McKeown, “Multiblock PLS model for group corticomuscular activity analysis in Parkinson disease,” in Conference Record of the Forty Fourth Asilomar Conference on Signals, Systems and Computers (ASILOMAR’10), Nov. 2010, pp. 1375-1379. In both Chapters 6 and 7, I was responsible for technical literature survey, development of analysis techniques, numer- ical simulation, and analysis and compilation of results. The EEG data analyzed in these two chapters were collected by Samantha Palmer. I prepared the manuscripts with comments and suggestions from Dr. Z. Jane Wang and Dr. Martin J. McKeown. The neurological interpretation of EEG results was written with the guidance of Dr. Martin J. McKeown. Other publications resulting from this PhD work are: • G. Tropini, J. Chiang, Z. J. Wang, and M. J. McKeown, “Partial directed coherence- based information flow in Parkinson’s disease patients performing a visually-guided motor task,” in Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC’09), Sep. 2009, pp. 1873-1878. • J. Chiang, Z. J. Wang, and M. J. McKeown, “EEG source extraction by autoregressive source separation reveals abnormal synchronization in Parkinson’s disease,” in Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC’09), Sep. 2009, pp. 1868-1872. • J. Chiang, G. Tropini, Z. J. Wang and M. J. Mckeown, “Broad-Band directional EEG v connectivity changes during motor preparation predicts Parkinson’s disease severity”, in The XIXWFNWorld Congress on Parkinson’s Disease and Related Disorders, Shang- hai, China, 2011. • G. Tropini, J. Chiang, Z. J. Wang and M. J. Mckeown, “Broad-Band directional EEG connectivity changes during motor preparation predicts Parkinson’s disease severity”, in The XIXWFNWorld Congress on Parkinson’s Disease and Related Disorders, Shang- hai, China, 2011. All experimental procedures in this work were reviewed and approved by the University of British Columbia Clinical Research Ethics Board under Certificate Numbers H04-70291 and H05-70555. vi Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii List of Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Brain Connectivity Analysis with EEG and Challenges . . . . . . . . . . . . 3 1.2 Related Works on EEG and EMG Connectivity Analysis Techniques . . . . . 4 1.2.1 Brain Connectivity Measures . . . . . . . . . . . . . . . . . . . . . . 4 1.2.2 Brain Source Extraction . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2.3 EEG-EMG Coupling . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3 Assessment of Abnormal Connectivity in Parkinson’s Disease Subjects . . . . 9 1.3.1 Background on Parkinson’s Disease . . . . . . . . . . . . . . . . . . 9 1.3.2 Modulation of Connectivity by PD . . . . . . . . . . . . . . . . . . . 10 1.4 Motivation and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.5 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 vii 2 Sparse Partial Directed Coherence (PDC) for EEG Analysis . . . . . . . . . . . 16 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.1 Sparse mAR Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.2 Partial Directed Coherence . . . . . . . . . . . . . . . . . . . . . . . 22 2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.3.1 Simulated Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.3.2 Case Study: EEG Analysis . . . . . . . . . . . . . . . . . . . . . . . 25 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3 Altered Directional Connectivity in Parkinson’s Disease Using Sparse PDC . . 30 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2.2 Sparse mAR-PDC Connectivity Analysis . . . . . . . . . . . . . . . 33 3.2.3 Spectral Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.3 Experimental Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.3.1 Subjects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.3.2 Behavioral Paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3.3 Data Collection and Analysis . . . . . . . . . . . . . . . . . . . . . . 36 3.3.4 Movement Preparation (PRE) and Movement (MOV) Phases . . . . . 38 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.4.1 Behavioral Measurements . . . . . . . . . . . . . . . . . . . . . . . 40 3.4.2 Frequency Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.4.3 Connectivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4 Prediction of PD Severity Using Sparse PDC Connectivity Strengths . . . . . . 58 viii 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.2.1 Feature Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.2.2 Cross Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.3 Experimental Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5 Sparse Group mAR Model for Assessing Effects of Galvanic Vestibular Stim- ulation on Brain Connectivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.2.1 Spectral Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.2.2 Sparse Group mAR Model . . . . . . . . . . . . . . . . . . . . . . . 73 5.3 Experimental Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.3.1 Subjects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.3.2 Stimulus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.3.3 Primary Study Protocol . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.3.4 Data Collection and Preprocessing . . . . . . . . . . . . . . . . . . . 76 5.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.4.1 Threshold Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.4.2 EEG Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 6 Generalized Multivariate Autoregressive Framework for EEG Source Separa- tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 6.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 6.2.1 EEG Signal Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 ix 6.2.2 Maximum Likelihood Estimation . . . . . . . . . . . . . . . . . . . . 89 6.2.3 Group Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 6.2.4 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.3 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.4 Application to Real EEG Data - A Case Study . . . . . . . . . . . . . . . . . 108 6.4.1 Experimental Procedures and Data Acquisition . . . . . . . . . . . . 108 6.4.2 Data Preprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 6.4.3 EEG Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 6.5 Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 7 A Multiblock PLS Model for Corticomuscular Interactions . . . . . . . . . . . 122 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 7.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 7.2.1 PLS Model and Estimation . . . . . . . . . . . . . . . . . . . . . . . 125 7.2.2 Application of mbPLS and EEG/EMG Feature Generation . . . . . . 130 7.2.3 Assessment of Significance . . . . . . . . . . . . . . . . . . . . . . . 134 7.3 Experimental Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 7.3.1 EEG Data Acquisition and Preprocessing . . . . . . . . . . . . . . . 136 7.3.2 EMG Data Acquisition and Preprocessing . . . . . . . . . . . . . . . 137 7.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 7.4.1 EEG Spectrogram Components Corresponding to EMG and Behav- ioral Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 7.4.2 PDC Components Corresponding to EMG and Behavioral Data . . . . 141 7.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 7.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 8.1 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 x 8.2 Future Research Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 A Estimation of Multiblock PLS . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 xi List of Tables Table 3.1 Summary of patient details. . . . . . . . . . . . . . . . . . . . . . . . . . 35 xii List of Figures Figure 2.1 Examples of penalty functions. . . . . . . . . . . . . . . . . . . . . . . . 21 Figure 2.2 PDC networks of simulated data. . . . . . . . . . . . . . . . . . . . . . . 24 Figure 2.3 PDC of selected connections. . . . . . . . . . . . . . . . . . . . . . . . . 26 Figure 2.4 Head plots showing significant connections in Gamma band (36-60Hz). . 28 Figure 2.5 Average4PDCregular and4PDCsparse for the connection from CZ to F8. 29 Figure 3.1 Joystick task schematics with non-predictive cue. . . . . . . . . . . . . . 36 Figure 3.2 Headplot of electrode regions. Electrode placement is based on the 10-20 system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 Figure 3.3 Spectrograms for all groups. . . . . . . . . . . . . . . . . . . . . . . . . 42 Figure 3.4 Spectra during PRE and MOV phases for the five regions. . . . . . . . . . 44 Figure 3.5 PDCograms of FCentral↔ Central and LSm↔ Occipital for all groups. . 46 Figure 3.6 PDCograms of RSm↔ Occipital and Central↔ Occipital for all groups. 47 Figure 3.7 PDC Spectra for the PRE phase of the task. . . . . . . . . . . . . . . . . 48 Figure 3.8 PDC Spectra for the MOV phase of the task. . . . . . . . . . . . . . . . . 49 Figure 3.9 Correlation between low-gamma band PDC in the Central → FCentral connection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 Figure 4.1 Illustration of cross validation procedure. . . . . . . . . . . . . . . . . . 62 Figure 4.2 Joystick task schematics with predictive warning cue. . . . . . . . . . . . 63 Figure 4.3 Prediction performance for the movement preparation phase. . . . . . . . 65 Figure 4.4 Prediction performance for the movement execution phase. . . . . . . . . 66 xiii Figure 4.5 Features selected by LASSO for the MOV-PREP phase and their respec- tive number of times being selected in the jackknife procedure. . . . . . . 67 Figure 4.6 Features selected by LASSO for the MOV-EXEC phase and their respec- tive number of times being selected in the jackknife procedure. . . . . . . 68 Figure 4.7 Predicted UPDRS values of each subject group using connectivity from the MOV-PREP phase. . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 Figure 5.1 Characteristics of stimuli delivered. . . . . . . . . . . . . . . . . . . . . 76 Figure 5.2 Power spectra during application of GVS. . . . . . . . . . . . . . . . . . 79 Figure 5.3 PDC spectra during application of GVS. . . . . . . . . . . . . . . . . . . 80 Figure 6.1 Estimation procedure for group connectivity analysis. . . . . . . . . . . . 92 Figure 6.2 Spatial characteristics of a source dipole. . . . . . . . . . . . . . . . . . . 94 Figure 6.3 PDC estimation results from an example simulation run. . . . . . . . . . 99 Figure 6.4 Scalp projection maps of the estimated sources given by GmAR-SS . . . 100 Figure 6.5 Effects of sensor noise on mixing matrix and PDC estimation. . . . . . . 103 Figure 6.6 Effects of biological noise on mixing matrix and PDC estimation. . . . . 104 Figure 6.7 Effects of sample length on mixing matrix and PDC estimation. . . . . . 105 Figure 6.8 Effects of Gaussianity of source innovations on mixing matrix and PDC estimation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 Figure 6.9 Scalp projection maps of the 6 source dipoles recovered by GmAR-SS . . 107 Figure 6.10 PDC of the selected estimated sources and simulated sources. . . . . . . . 109 Figure 6.11 PDC of all the 6 estimated sources given by GmAR-SS and the signifi- cance threshold determined using surrogate data. . . . . . . . . . . . . . 110 Figure 6.12 Connectivity patterns of simulated and GmAR-SS estimated sources . . . 111 Figure 6.13 Squeezing task. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 Figure 6.14 Groupings of electrodes in six brain regions . . . . . . . . . . . . . . . . 113 Figure 6.15 Dendrogram and similarity matrix. . . . . . . . . . . . . . . . . . . . . . 114 xiv Figure 6.16 Scalp projection maps of the estimated sources given by GmAR-SS from a typical PD subject. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 Figure 6.17 PDC connectivity patterns between EEG sources derived from spatial fil- tering. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 Figure 6.18 Average spectrograms of the EEG sources of normal subjects. . . . . . . 117 Figure 6.19 Average spectrograms of the EEG sources of PD subjects . . . . . . . . . 118 Figure 7.1 Multiblock PLS model. . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 Figure 7.2 EEG and EMG features. . . . . . . . . . . . . . . . . . . . . . . . . . . 132 Figure 7.3 Coherence between the raw EEG signal of the 5 region and the EMG signal of a typical normal subject. . . . . . . . . . . . . . . . . . . . . . 138 Figure 7.4 Spectrograms of raw EEG signals. . . . . . . . . . . . . . . . . . . . . . 138 Figure 7.5 Component 1 of PLS decomposition with EEG Spectrograms as predic- tors and EMG amplitude as responses. . . . . . . . . . . . . . . . . . . . 139 Figure 7.6 Component 7 of PLS decomposition with EEG Spectrograms as predic- tors and EMG amplitude as responses. . . . . . . . . . . . . . . . . . . . 140 Figure 7.7 Component 1 of mbPLS decomposition with EEG PDC as predictors and EMG amplitude as responses. . . . . . . . . . . . . . . . . . . . . . . . 143 Figure 7.8 Component 2 of mbPLS decomposition with EEG PDC as predictors and EMG amplitude as responses. . . . . . . . . . . . . . . . . . . . . . . . 144 Figure 7.9 Component 8 of mbPLS decomposition with EEG PDC as predictors and EMG amplitude as responses. . . . . . . . . . . . . . . . . . . . . . . . 145 xv List of Acronyms AIC Akaike information criterion ANOVA analysis of variance AR autoregressive BG basal ganglia BIC Bayesian information criterion BSS blind source separation DBS deep brain stimulation DTF directed transfer function EEG electroencephalogram EMG electromyogram ERS event-related synchronization ERD event-related desynchronization FDR false discovery rate fMRI functional magnetic resonance imaging GC Granger causality GPi globus pallidus interna GVS galvanic vestibular stimulation ICA independent component analysis LARS least angle regression LFP local field potential LQA local quadratic approximation LS least square xvi mAR multivariate autoregressive MEG magnetoencephalogram ML maximum likelihood MSC magnitude-squared coherence PCA principal component analysis PD Parkinson’s disease PDC partial directed coherence PET positron emission tomography PLS partial least square ROI region of interest SL synchronization likelihood SBNR signal-to-biological noise ratio SNR signal-to-sensor noise ratio STN subthalamic nucleus TMS transcranial magnetic stimulation UPDRS Unified Parkinson’s Disease Rating Scale VA ventroanterior VL ventrolateral xvii Acknowledgments I would like to take this opportunity to convey my great appreciation to my supervisors, Dr. Z. Jane Wang and Dr. Martin J. McKeown, whose patient guidance, persistent encouragement and deep insight in the research area have helped me immeasurably throughout the course of my thesis research. This thesis would never have been written without their assistance. I would like to thank my thesis examination committee members for their time and ef- forts. I would also like to acknowledge Giorgia Tropini and Dr. Edna Ty from UBC Brain Research Centre for providing us with the data and helpful discussions regarding the medical interpretation of the results. I am deeply indebted to my family and C.H. for their constant support and immense encouragement over the years. This work was supported by the Canadian Natural Sciences and Engineering Research Council (NSERC) Strategic Project Grant (STPGP 365208-08) and Postgraduate Scholarship (PGS 362898-08) and by the Canadian Institutes of Health Research (CIHR) Grant. xviii Chapter 1 Introduction A fundamental question in neuroscience research is how the brain is able to perceive, process and act upon the tremendous amount of information flowing across in a seemingly effortless manner. While the anatomical structure of the brain is relatively well understood, the under- lying mechanism that coordinates the various components to act together synergistically is yet to be elucidated. However, the problem of trying to understand the brain and establish the relationships between individual components is highly complex as it involves billions of neurons interlinked through trillions of connections. Research efforts in this area have looked at different levels of the brain, ranging from the molecular level which includes identify- ing anatomical connections between neurons and mechanisms for relaying information be- tween neurons to the system level which models the flows of information between neuronal populations and formation of cognitive functions. A major part of existing works focus on identifying the functional specialization of a particular brain region or functional segregation based on the physiological responses caused by external stimuli or localized cortical legions. However, most cognitive tasks, as simple as grasping a ball or pushing a button, involve the coherent activation of several functionally specialized regions. Characterizing brain activity in terms of functional specialization thus only provides a limited view of brain function and does not explain how different regions communicate with each other. With the availability of multichannel neural signals, researchers start to analyze these 1 data in terms of the long-range interactivity between different functionally specialized neu- ronal populations or functional integration. This has led to the emergence of the concept of brain connectivity. The study of brain connectivity has provided not only a system-level view of brain functioning in the normal state but also an explanatory framework for pathological conditions such as autism [1], Alzheimer [2] and Parkinson’s disease [3, 4], which are often manifested in the form of disrupted or altered connectivity patterns. Brain connectivity can be characterized in two different forms, functional and effective connectivity. As defined in [5], functional connectivity is the temporal correlation between two spatially remote neuro- physiological events, which is typically expressed as deviation from statistical independence. Effective connectivity refers to the influence that one neural entity exerts over another directly or indirectly. The inference of effective connectivity involves establishing causal relationships between different neural entities via mathematical modeling of their neural activity. In this thesis, we are particularly interested in identifying effective connectivity as it pro- vides an explanation for the mechanism underlying brain function, complementing the con- ventional functional connectivity analysis which is only able to identify correlative patterns of interactions between neural entities. Using modeling and computational theories from signal processing, novel connectivity analysis techniques will be developed and applied to normal and Parkinson’s disease (PD) subjects during dynamic movements in an attempt to assess the effects of disease and medication on motor-related brain connectivity. We further extend the connectivity analysis from the cortico-cortical level to cortico-muscular level by examining the relationship between simultaneously recorded cortical and muscular activity. The ultimate goal of this research work is to provide a set of computational tools which not only enhance our understanding of how the brain works in the normal state but also provide an explanation for what makes it go astray and potential solutions to make it function correctly again. The rest of this chapter is organized as follows: Section 1.1 outlines the advantages and challenges associated with measuring brain connectivity using scalp electroencephalogram (EEG). Section 1.2 reviews the most commonly used methods for measuring EEG brain connectivity. Section 1.3 provides an introduction to the PD and its effects on brain activity. 2 The research objectives are presented in Section 1.4. Finally, Section 1.5 presents the thesis overview. 1.1 Brain Connectivity Analysis with EEG and Challenges Various techniques and conceptual frameworks from physics, mathematics and engineer- ing have been explored to measure brain connectivity. The choice of analysis technique depends on the statistical characteristics of the data being analyzed including spatial and temporal resolution and level of description, e.g., activity of individual neurons, activity of neuronal populations or large-scale activations [6]. In particular, noninvasive neuroimaging techniques, including functional magnetic resonance imaging (fMRI), positron emission to- mography (PET), magnetoencephalogram (MEG) and EEG, are becoming prevalent in recent neuroscience studies as they allow assessment of large-scale cortical interactions in vivo while subjects perform cognitive/motor tasks. Among these, scalp EEG is one of the most widely used noninvasive techniques. Scalp EEG measures brain activity by placing multiple elec- trodes on the scalp surface. The measured signals reflect the summation of the synchronous activity of a large number of spatially aligned neurons (thousands or millions) in the brain. EEG is more sensitive to the activity generated in the superficial layers of the cortex and activity of deeper brain structures is more difficult to detect [7]. Scalp EEG offers several advantages over other recording techniques. It is noninvasive and is relatively more tolerant of subject movement. More importantly, EEG also has high temporal resolution (on the level of milliseconds) compared of fMRI and PET which have time resolution between seconds and minutes, thus allowing the detection of transient ac- tivity [7]. Nonetheless, the assessment of brain connectivity using EEG recordings is con- founded by two main factors. One is the poor signal-to-noise ratio of EEG, where the ‘noise’ includes artifacts associated with eye movements, blinks and muscle noise, and non-task re- lated background EEG activities. While there are several effective algorithms to remove eye and muscle artifacts, background EEG activities are much more difficult to characterize and isolate, which makes the extraction of task-related brain activities challenging. Another ma- 3 jor issue plaguing the analysis of EEG signals is the problem of volume conduction [6]. It arises from the fact that brain tissues act as a conductor and as an electric field gets generated by the firing neurons (or dipoles), it gets propagated radially to the scalp and picked up by multiple measurement sensors. In other words, the signal measured from an EEG electrode does not exclusively represent the activity of one local neural source, but rather the superpo- sition of several active sources throughout the brain [8]. The poor spatial resolution of scalp EEG complicates the correct interpretation of connectivity results. 1.2 Related Works on EEG and EMG Connectivity Analysis Techniques In this section, we will review the most common techniques that have been used for analyzing EEG and electromyogram (EMG) connectivity. 1.2.1 Brain Connectivity Measures To provide a quantitative measure of connectivity between oscillatory signals, several time series analysis techniques have been explored in the literature. They vary in terms of direct- edness (symmetric vs. directional), character of interdependence (linear vs. nonlinear) and the number of signals being considered (bivariate vs. multivariate). Here, we will restrict ourselves to highlight the ones that are more widely used in EEG connectivity studies. More detailed surveys can be found in [9] and [10]. The most classical connectivity measure is spectral coherence which measures the linear correlations between two signals as a function of frequency. Spectral coherence has been used extensively in EEG studies to investigate cortical synchrony during different cognitive tasks and disrupted brain connectivity in pathological conditions. Despite its popularity, the coherence technique has several limitations. First, it is a bivariate technique, meaning that only two signals are considered at a time. By doing pairwise analysis, coherence technique ignores possible critical influences by other neural sources and is more sensitive to the pres- ence of volume conduction. Moreover, it assumes the symmetry of the connection and is 4 unable to provide directionality of information flow between neural sources. Techniques that measure statistical interdependence based on information theory have also been explored. The most common used one is mutual information, which is given by MIXY = ∫ ∫ pXY (x,y) log pXY (x,y) pX(x)pY (y) dxdy, (1.1) where pXY is the joint probability density function of X and Y and pX and pY are the marginal probability density functions of X andY , respectively [11]. It essentially tells how much extra information one gets from one signal by knowing the outcomes of the other one. It is a time domain measure which allows the quantification of nonlinear dependencies between two neu- ral sources. Despite the simplicity of its definition, the estimation of mutual information from real data is not straightforward. The most common way to estimate the probability density functions from observed signals is to obtain the histograms of the series of outcomes. How- ever, to get an accurate estimate of the mutual information using histogram-derived probabili- ties, a large number of data samples is required, which is often prohibitive in real applications. In addition, similar to the coherence method, mutual information is a symmetric measure and thus does not provide information on the directionality of interactions [9]. To deal with the limitations of the aforementioned techniques, several directional interac- tion analysis techniques have been proposed. One technique that has been commonly used to evaluate causal relations is Granger causality (GC). The basic concept of GC is that if the prediction of a time series, x(t), can be improved by incorporating the knowledge of a second time series, y(t), then y(t) is said to have a causal influence on x(t). This idea was originally proposed by Norbert Wiener and later formalized by Clive Granger in the context of a linear regression framework, where the causality is assessed by comparing the fitting of a univariate autoregressive (AR) model to the signals with that of a bivariate model [9]. Extending the concept of GC, several multivariate, frequency-domain causality measures are proposed, including directed transfer function (DTF) [12] and partial directed coherence (PDC) [13]. These techniques involve modeling multichannel signals by a single multivariate 5 autoregressive (mAR) process, which is formulated as follows: x(t) = P ∑ p=1 Apx(t− p)+e(t), (1.2) where x(t) is anM-dimensional observation vector at time t and e(t) represents anM-dimensional white Gaussian noise. The matrix Ap denotes the mAR coefficient matrix, where the element Ap(i, j) describes the causal influence of signal j on signal i after p time points. This time- domain representation can be translated to frequency domain by taking Fourier transform of the mAR coefficients. The transformed mAR coefficients provide a quantification of the causal influences as a functional of frequency. Compared to conventional bivariate tech- niques, DTF and PDC give more accurate estimations of causality for multivariate time series as they account for all available correlation information from the entire data set and are more robust to the effects of volume conduction compared to the conventional coherence method. Furthermore, the spectral information provided by DTF and PDC is particularly useful for EEG studies since EEG signals often exhibit frequency-specific oscillatory activities. The difference between PDC and DTF is subtle; PDC is normalized with respect to all the out- flows from a given source channel whereas DTF is normalized with respect to all the inflows to a given destination channel. In [13], it was further shown that PDC is able to break up the interactions between channels into their direct pair-wise components thus allowing eas- ier interpretation of results [13, 14]. Examples of applications of PDC to biomedical signals include the study of brain networks in rats during different behavioral states [15] and the iden- tification of oscillatory brain interactions in human during object recognition [16]. Although GC-based techniques can only capture linear interdependencies, they are computationally efficient and offer a straightforward characterization of effective connectivity [17, 18]. How- ever, in order to obtain accurate estimate of the connectivity pattern, parameters such as mAR model order, number of channels, number of data samples should be carefully chosen [14, 19, 20]. 6 1.2.2 Brain Source Extraction One major problem plaguing EEG analysis is volume conduction, which affects the spatial resolution of recorded brain activity. To address this issue, several studies have proposed studying the brain connectivity in the “source space” which involves estimating the underly- ing neural sources from scalp EEG signals. Determining neural sources from EEG signals is a mathematically ill-conditioned inverse problem, and it has no unique solution without prior knowledge or strong statistical assumptions. One common approach to solving the source identification problem is the dipole modeling technique [21], which aims to find cur- rent dipoles that best describe the observed EEG signals. The dipole localization process involves iterative estimation of 3 location parameters, 2 orientation parameters and one mo- ment (strength) parameter for each dipole, until the scalp potential predicted at the measure- ments sites based on the assumed head model match with the actual EEG recordings. To limit the number of possible solutions, a priori information based on multimodal integration tech- niques such as co-registration of EEG with MRI, fMRI or PET, is often incorporated into the source model. However, the successful estimation of sources depends on many factors, in- cluding models of the generators of the electrical activity, knowledge of the number of active dipoles, and the shape and conductivity of the head model [22]. Another possible approach to finding brain sources is to use blind source separation (BSS) algorithms. BSS aims to extract a set of underlying sources from observed signals by mak- ing certain assumptions about the statistical properties of the source signals. One of the most widely used BSS algorithms in neuroscience studies is independent component analy- sis (ICA) . Fundamentally, ICA assumes the sources are mutually independent. Most ICA algorithms further assume that observed signals are instantaneous linear mixtures of the un- derlying sources and their temporal ordering is irrelevant [23]. ICA is often used to identify and remove eye movement and muscle noise artifacts in EEG signals, as these artifacts are typically independent from the brain activities of interest. However, in the context of con- nectivity analysis, as we try to determine the causal relationships between EEG sources, it 7 is implicitly assumed that the sources are not only temporally correlated within themselves but also cross-correlated between each other. This contradicts the fundamental assumption of ICA. 1.2.3 EEG-EMG Coupling Human movement control involves complex interplay between the brain and the musculature. To explore the functional role of cortical activity in the human motor system, a possible approach is to assess the interactivity between cortical activity and simultaneously record muscular activity or cortico-muscular coupling during muscle contraction [24, 25]. Previous studies using coherence as a measure of cortico-muscular coupling show beta-range (15-40 Hz) synchronization between EMG of the contracting muscles and MEG [26, 27] or EEG [28] of the primary sensorimotor cortex during sustained isometric (static) contractions and lower- frequency (2-14 Hz) synchronization during phasic voluntary movement [29]. Pathological alternations in cortico-muscular coherence are observed in patients exhibiting motor-related disorders [30, 31]. Despite the common use of the coherence method for assessing EEG-EMG coupling, it suffers from several limitations. First, the common formulation of coherence relies on pairwise comparisons yet it has been suggested that the mapping between the brain and mus- culature during movements is many-to-many, as opposed to one-to-one [32]. Hence, the utilization of multivariate techniques to look at several EEG and EMG channels simultane- ously may enhance the sensitivity in detecting cortico-muscular coupling. Second, an EEG recording represents the superposition of several ongoing brain activities and are often con- taminated by eye and movement artifacts and measurement noise. Only a small fraction of the EEG signals is related to the control mechanism underlying the performed motor task [33]. Hence, the conventional approach of computing coherence between raw EEG and EMG typically yields a very low coherence value and sometimes not being distinguishable at all for some percentage of population. Nonetheless, based on the recent findings of EEG source identification, one could expect that with a proper transformation of the raw EEG recordings, 8 a waveform corresponding to the actual drive to the muscle can be recovered and presumably yielding a stronger correlation with the EMG signals of active muscles. 1.3 Assessment of Abnormal Connectivity in Parkinson’s Disease Subjects 1.3.1 Background on Parkinson’s Disease Parkinson’s disease is a neurodegenerative disorder resulted from the death of cells in the substantia nigra that are responsible for producing dopamine - a neurotransmitter chemical that acts as a messenger between the cells of the brain that control movements [34]. PD is characterized by motor symptoms such as tremor at rest, rigidity, slowness of movement (bradykinesia) and difficulty with walking and gait. As the disease progresses, non-motor symptoms such as cognitive and behavioral alterations including dementia may also appear. PD is the second most common neurodegenerative disorder after Alzheimer’s disease and it is estimated that more than 100,000 Canadians are diagnosed with PD [35]. The severity of PD is assessed clinically using a widely used metric, the Unified Parkin- son’s Disease Rating Scale (UPDRS). The scale consists of 4 components assessing behavior and mood, activities of daily living, motor symptoms and complication of therapy. These first 3 components contain 44 questions, each measured on a 5-point scale (0-4). The ratings of these three components are combined to give a total UPDRS score with 0 being not affected and 176 being most severely affected. In adjunct therapy, the fourth component containing 11 questions with a scale ranging from 0 to 23 is included [36]. The UPDRS is often used in con- junction with other measures (e.g., Hoehn and Yahr score) to evaluate motor symptoms and clinical status in PD patients. Despite the wide utilization of UPDRS, the assessment is time- consuming and can only be performed by trained medical personnel. Furthermore, adjusting the dosage and monitoring the effects of PD medication may require frequent re-assessment. While there is currently no cure for PD, many of the motor symptoms can be controlled by pharmacological and surgical treatments [34]. The main families of drugs used to treat PD, 9 including levodopa and dopamine agonists, work by either temporarily replenishing the lost dopamine or imitating the action of dopamine in the brain. Surgical treatment via deep brain stimulation (DBS) is also possible [37]. DBS is particularly effective in treating patients who do not respond to medications or experience intolerable side effects of their medications. DBS operation involves implantation of a medical device called brain pacemaker, which sends high frequency electrical pulses to specific areas of the brain critical for functional control of motor activity. While these two approaches have been shown effective in suppressing the symptoms, they can result in serious side effects such as compulsive behaviors, cognitive decline, depression and psychosis [34]. 1.3.2 Modulation of Connectivity by PD The Parkinsonian state is also manifested in the form of abnormal brain activity. The most prominent characteristics in brain activity of PD patients is the increase in oscillations in the beta band (13-30 Hz) [38, 39]. Several studies report that local field potential (LFP) recorded from PD patients undergoing a DBS procedure consistently show increased beta band power at various sites within the basal ganglia loop including subthalamic nucleus (STN) and globus pallidus interna (GPi). Beta synchrony is normally associated with the movement-related activity where it is more prominent during holding contractions and is attenuated shortly before and during voluntary movements. This has lead to the hypothesis that the excessive beta band synchronization may be antikinetic in nature and is causally related to the difficulty in initiating movements seen in PD patients [40]. The administration of levodopa medications and DBS are shown effective in suppressing the abnormalities in beta band synchronization in PD patients [41] and the degree of levodopa-induced beta band power suppression correlates positively with clinical improvement in motor performance. Abnormal oscillatory activity in PD is not limited to LFP of STN and GPi but is also ob- served at the cortical level via noninvasive neuroimaging recordings. A number of previous studies have investigated the relationship between the activity of STN and cortex via concur- rent LFP and scalp EEG recordings. Most of these studies utilize coherence as a measure of 10 correlation and only a limited number of cortical regions are measured. Later studies start to look at altered cortical activity in PD subjects measured by EEG and MEG. In [42], EEG recordings were taken in PD patients after overnight withdrawal from medication and again after 1 hour after levodopa administration. Levodopa-related attenuation of the alpha (8-12 Hz) and beta band power over the cortical motor area is observed and is correlated with im- provement in speed of movement. Other studies using EEG recordings have reported similar observation. These result are in line with observations made using LFP recordings of STN, suggesting that EEG is suitable for studying PD. EEG has several advantages over the conventional LFP technique for studying brain ac- tivity in PD. In majority of the studies on abnormal brain activity in PD, LFP recordings are only available from PD patients undergoing the DBS surgery and this precludes the healthy subjects and PD patients at less advanced stages of the disease from the study. Given the noninvasive nature of EEG, EEG measurements can thus be taken from both PD and healthy subjects, enabling the comparison of brain activity between these two subject groups. Further- more, EEG provides a wider coverage of the brain and thus is useful for studying long-range cortical interactions. A limited number of EEG studies have investigated cortico-cortical functional connectivity in PD patients, mainly using coherence as a measure of coupling be- tween cortical regions. A more complete analysis of pathological effects of PD on directional connectivity patterns and comparison with healthy subjects is yet to be done. 1.4 Motivation and Objectives The concept of brain connectivity has provided a new perspective to the understanding of the mechanism underlying brain function in both normal and disease states, complementing the conventional approach of focusing on region-specific changes in amplitude or spectral distribution. The wide use of noninvasive recording techniques has further accelerated this paradigm shift in ways of analyzing neurophysiological data. In particular, scalp EEG has been used in a growing number of studies to study brain connectivity in pathological condi- tions. Recent EEG studies on PD reveal the association between motor deficits observed in 11 patients withdrawn from medication and alterations in connectivity patterns. Most of these studies are based on coherence as the measure of functional connectivity and fail to pro- vide information about the causal relationships between the involved cortical regions. In this thesis, we will expand on the previous studies by using mAR model-based connectivity mea- sures to infer causal dependencies as the mAR model allows the simultaneous modeling of all EEG channels, giving more reliable estimations of causality than conventional bivariate techniques. However, the correct inference of brain connectivity from EEG recordings using mAR- based measures is complicated by several factors. In addition to the poor signal-to-noise ratio and effects of volume conduction inherent in EEG signals as described in Section 1.1, another challenge comes from statistical instability during the estimation of mAR model parameters. In many real EEG applications, the number of acquired EEG channels can vary from 20 to up to 128 channels. Furthermore, in order to achieve the desired frequency resolution of the connectivity measure, a high mAR model order is typically required. These two factors result in a substantial increase in the number of mAR parameters to be estimated, which in turn affects the statistical stability of ordinary maximum likelihood estimators when only limited data samples are available. Lastly, when applying connectivity analysis to multiple subjects for group analyses, the presence of inter-subject variability is often overlooked or not dealt with explicitly. Conventional approaches to handling inter-subject variability involve either averaging/concatenating the data across subjects and performing connectivity analysis to the averaged/concatenated data assuming that the connectivity pattern is identical across subjects, or performing connectivity analysis on each subject separately and later integrating the results across subjects using statistical tests. These two approaches represent two extremes in finding group connectivity networks: the former allows the identification of common features across subjects while the latter preserves the subject specificity which may be used to correlate with clinical or demographic measures (e.g., UPDRS or task performance) of the subjects. Further exploration in group analysis techniques to achieve a balance between these two extremes is warranted. 12 The goal of this research work is to develop novel signal processing methods to address the aforementioned technical challenges arising when inferring brain connectivity from EEG recordings. Specifically, the main technical objectives are: 1. Develop a new computational model for assessing EEG connectivity which incorpo- rates the concept of sparsity to circumvent the statistical instability in mAR model estimation. 2. Develop blind source separation methods which simultaneously model the underlying sources and connectivity between them in the presence of volume conduction. 3. Address the issue of inter-subject variability in group-level connectivity analysis and source modeling. 4. Propose a multi-subject, multi-modal data analysis framework for investigating rela- tionship between cortical and muscular activity in PD patients during dynamic motor tasks. The developed connectivity analysis techniques are applied to EEG data collected from healthy subjects and PD patients, including ‘on’ medication and ‘off’ medication states dur- ing the performance of various motor tasks. Furthermore, a novel algorithm which extracts informative features from EEG connectivity measures to predict PD severity is proposed. While UPDRS has been the gold standard assessment tool for characterizing impairments in PD patients, it is time-consuming and can only be performed by trained medical personnel. Previous studies have demonstrated the possibility of using features from speech recordings to predict PD severity [43, 44]. Our approach of using brain connectivity as predictive fea- tures of disease severity can have useful clinical implications in understanding the effects of disease progression on connectivity modulation. 1.5 Thesis Overview In Chapter 2, we propose a sparse mAR-based PDC technique where PDC estimates are de- rived from mAR models with sparsity constraints imposed on subject-level mAR coefficients. 13 The sparse mAR parameters are estimated using penalized regression which automatically carries out variable selection, in effect reducing the number of free parameters during estima- tion. The proposed technique is applied to both simulated data and real EEG recordings to demonstrate the improved estimation performance of the proposed technique relative to the conventional non-sparse approach. In Chapter 3, the sparse mAR-based PDC technique developed in Chapter 2 is applied to PD and normal subjects performing a visually guided motor task. The connectivity results are compared with the results obtained utilizing traditional Fourier analysis to evaluate task- dependent frequency spectra modulation. In Chapter 4, we investigated the relationship of the functional connectivity during move- ment preparation and execution with clinical measures of motor function in PD subjects. Based on the connectivity strengths measured using sparse mAR-based PDC, we propose a novel feature selection technique to identify brain connections whose connectivity strengths are able to both differentiate PD subjects from normal subjects and accurately predict UPDRS scores. In Chapter 5, we extended the sparse mAR model proposed in Chapter 2 and developed a group sparse mAR model which incorporate the Group LASSO (gLASSO) method into the estimation of mAR coefficients. The gLASSO enforces the mAR models of the subjects within a group to be structurally identical across subjects, but the connection coefficients are allowed to vary from subject to subject. The group sparse mAR model is applied to EEG data collected from healthy subjects undergoing the galvanic vestibular stimulation. In Chapter 6, we proposed a state-space framework to jointly model the instantaneous mixing effects of volume conduction and the causal relationships between underlying brain sources. We modeled the source activity by a generalized mAR process with possibly non- Gaussian noise. A maximum likelihood estimation approach is developed which allows si- multaneous estimation of both the mixing matrix and AR model parameters directly from scalp EEG. The proposed technique was verified with simulated EEG data generated using the single-shell spherical head model and compared to conventional two-stage connectivity 14 estimation approaches. In addition, the proposed technique was applied to EEG data collected from normal and PD subjects performing a right-handed force-tracking task with differing amounts of visual feedback. In Chapter 7, we proposed a multiblock partial least square (mbPLS) framework for mod- eling the correlative relationship between concurrent EEG and EMG signals at a group level. The mbPLS framework incorporates a two-level hierarchical structure into the ordinary two- block partial least square (PLS) with subject-specific details separated into individual data blocks at the sub-level and the common group patterns at the super-level. We further ex- tended mbPLS to enable modeling of 3-dimensional matrices: time-frequency-EEG channel and time-frequency-connection utilizing PDC. The proposed 3-dimensional mbPLS frame- work was applied to concurrent EEG and EMG data collected from ten normal subjects and nine patients with mild-moderate PD performing a dynamic squeezing task. Finally, the main contributions of the thesis and suggestions for future research are pre- sented in Chapter 8. 15 Chapter 2 Sparse Partial Directed Coherence (PDC) for EEG Analysis 1 2.1 Introduction In neurobiology, there has been increasing interest in identifying functional connectivity between brain regions. Such connectivity is believed to provide an integrating framework for a variety of complex brain functions. Several mathematical methods have been ex- plored in the literature to provide a quantitative measure of brain connectivity based on elec- troencephalogram (EEG) data, including correlation, coherence and Granger causality (GC) [9, 13]. Among them, spectral coherence, probably the most popular one, has been used extensively in EEG studies to investigate topics such as cortical synchrony during cognitive tasks and disrupted brain connectivity in pathological conditions. However, coherence tech- nique has several limitations. First, it is a bivariate technique, meaning that only two signals are considered at a time. In situations where large number of brain regions are simultaneously interacting with each other, pairwise analysis essentially ignores possible critical influences 1A version of this chapter has been published as: J. Chiang, Z. J. Wang, and M. J. McKeown, “Sparse multivariate autoregressive (mAR)-based partial directed coherence (PDC) for electroencephalogram (EEG) analysis,” in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP’09), Apr. 2009, pp. 457-460. c© 2009 IEEE. http://dx.doi.org/10.1109/ICASSP.2009.4959619 16 by other regions and thus may lead to misleading results. Moreover, coherence is unable to identify the direction of information flow between cortical regions [9]. Another technique that has received great attention recently is GC as it provides a measure of causal relation between two time series. Its basic idea is that if the prediction of one time series can be improved by incorporating the knowledge of a second one, then the second time series is said to have a causal influence on the first. Granger formalized this idea in the context of linear regression models where the causal influence is determined by comparing the variance of residuals from a scalar autoregressive (AR) application to one signal x(t), with that from a bivariate AR application to x(t) and a potentially driving signal y(t) [45]. However, similar to coherence technique, GC also only allows multiple pairwise analyses when considering the multichannel case [13]. To overcome these limitations, partial directed coherence (PDC) was later proposed for EEG connectivity studies [13]. PDC extends the concepts of coherence and GC to multi- channel time series and allows the simultaneous modeling of all channels with a multivariate autoregressive (mAR) model. This has the advantages of allowing differentiation of direct and indirect causal influences among interacting entities and giving more reliable estimations of causality than bivariate techniques. Moreover, since neural signals often exhibit frequency- specific oscillatory activity, the ability of providing spectral information of causal relations makes PDC an attractive tool for neuroscience studies. Examples of applications of PDC to real data include the study of brain networks in rats during different behavioral states [15] and the identification of oscillatory brain interactions in human during object recognition [16]. However, the computation of PDC poses several technical challenges in real applications, especially when the number of EEG channels of interest is relatively large in practice (e.g., 20). As noted in [9], the successful estimation of PDC depends on both the proper fitting of the mAR model to the data and also the accuracy of the specific parameter estimations, which in turn are dependent upon the number of channels, optimal mAR model order selection and sample sizes of training data [14]. In general, a higher model order allows more data dy- namics to be captured and gives a higher frequency resolution, but at the expense of a greater 17 number of parameters to be estimated. For instance, consider an EEG dataset consisting of 20 channels, each increment in the model order will result in an addition of 400 parameters to be determined. When only limited observations are available, such a large number of parameters may create statistical instability for ordinary maximum likelihood estimators [46]. Thus, it is of critical importance in practice to make a good trade-off between model complexity and estimation accuracy. In addition, the underlying assumption of full connections implicit in a regular mAR model used for PDC is questionable in EEG connectivity analysis, especially when the effects of volume conduction, which may erroneously suggest connections between channels, are account for. In other words, the connections between EEG nodes (or brain regions) may be considered a priori to form a sparse network. The above observations motivate us to incorporate a sparse mAR model into the current PDC technique. In an effort to resolve a similar concern in functional magnetic resonance imaging (fMRI) studies, Valdés-Sosa et al. [46] proposed modeling fMRI data using sparse mAR models in which the parameters are estimated using penalized regression. Penalized regression automatically shrinks small regression coefficients to zero during the estimation process. Thus, the resulting mAR coefficient matrix naturally have sparse structures, mean- ing that only small number of elements are non-zero. In this chapter, we propose a sparse mAR-based PDC technique in which PDC estimates are derived from sparse mAR coeffi- cient matrices given by penalized regression. Penalized regression effectively reduces the number of free parameters to be estimated, which is particularly important when the data is of limited sample size. We evaluate the performance of the proposed approach using both simulated data and real EEG data. This chapter is organized as follows: in Section 2.2, we first describe the regular and sparse mAR models and the parameter estimation techniques. The PDC function is defined in Section 2.2.2. In Section 2.3, the performance of regular mAR- and sparse mAR-based PDC are compared. Finally, we conclude in Section 2.4. 18 2.2 Methods We first review regular mAR models and least square-based parameter estimation techniques. We then introduce the concept of sparse mAR and present penalized regression methods for solving such sparse problems. Finally, the definition of PDC is given in Section 2.2.2. 2.2.1 Sparse mAR Model In a regular mAR model, the multivariate time series at each time point is represented as a linear, weighted sum of its previous time points, and it can be formulated as y(t) = P ∑ p=1 Apy(t− p)+e(t), (2.1) where the observation y(t) is a D-dimensional vector at time t, P denotes the order of the mAR model, and e(t) is a D-dimensional vector denoting multivariate white Gaussian noise. The mAR coefficient Ap is a D×D matrix, where element Ap(i, j) measures the influence node j exerts on node i after p time points. In the regression framework, Eqn. 2.1 can be rewritten as Z =XB+E (2.2) where Z = Y P+1:N = [ y(P+1), y(P+2), . . . , y(N) ]T ∈R(N−P)×D, X = [ Y P:N−1, Y P−1:N−2, . . . , Y 1:N−P ] ∈R(N−P)×DP, B = [ A1, A2, . . . , AP ]T ∈RDP×D, E = [ e(P+1), e(P+2), . . . , e(N) ]T ∈R(N−P)×D. Eqn. 2.2 can be solved using the maximum likelihood (ML) approach. Under the indepen- 19 dent and identically distributed (i.i.d.) white noise assumption ofE , the estimate of regression coefficients is obtained by minimizing the mean square error for each node i= 1,2, . . . ,D: β̂ i = argminβ i ‖(zi−Xβ i)‖2, (2.3) where zi and β i are the i-th column of Z and B, respectively. It is worth emphasizing that the performance of the ML estimator is highly dependent on the sample size N and the number of parameters to be estimated. In real applications, the available data points are often limited, which in turn leads to poor estimation accuracy. Furthermore, the estimated coefficients yielded by the least square (LS) approach in Eqn. 2.3 are typically non-zero, which makes neurobiological interpretation of results (e.g. identify- ing brain connectivity patterns in EEG studies) difficult. Such non-zero observation is also against the sparsity assumption in brain connectivity networks. To address these issues, a possible solution is to impose sparsity constraint on the mAR coefficients (i.e. Ap matrix) and perform variable selection using penalized regression meth- ods [46]. The basic idea of penalized regression is to maximize the likelihood while at the same time, penalize complex models. In terms of mathematical formulation, penalized re- gression can be expressed as the minimization of the penalized least square function: β̂ i = argminβ i ‖(zi−Xβ i)‖2+λ DP ∑ k=1 g(|βi(k)|), (2.4) where βi(k) is the k-th element of the regression coefficient vector β i and λ is the regular- ization parameter which controls the amount of penalization imposed on the solution. In the work, the value of λ is determined using Bayesian information criterion (BIC). BIC is widely used for model selection in linear regression and is defined as BIC=−2lnL+K ln(n), where L is the likelihood of the estimated model, K is the number of free parameters to be estimated and n is the number of observations. The model with the smallest value of BIC is selected. The function g(·) is the penalty function applied to each regression coefficient. Several dif- 20 −6 −4 −2 0 2 4 6 0 1 2 3 4 5 6 7 β g(β ) Ridge LASSO SCAD Figure 2.1: Examples of penalty functions. ferent penalty functions have been introduced in the literature, including ridge, LASSO and SCAD, and they are illustrated in Fig. 2.1. An overview of these penalty functions can be found in [47]. In this work, the popular LASSO penalization, g(|β j|) = |β j|, (2.5) is chosen due to its ability to automatically set small coefficients to zero, which effectively yields sparse solution. This special property, also known as the sparsity property, is particu- larly useful in variable selection problems [47]. To solve the optimization problem in Eqn. 2.4, the technique we use is the local quadratic approximation (LQA) algorithm, proposed by Fan and Li [47]. LQA first casts the prob- lem of penalized least square minimization presented in Eqn. 2.4 into a penalized likelihood maximization problem. It further addresses the issue of singularity at the origin that ex- ists in penalty functions such as LASSO and SCAD by locally approximating g(|β j|) with a quadratic function. The resulting penalized likelihood function becomes both differentiable and concave, and it can easily solved using gradient-based optimization algorithms such as Newton-Raphson. A detailed description of the LQA algorithm can be found in [47]. 21 2.2.2 Partial Directed Coherence PDC can be considered as the frequency-domain representation of GC. It involves the trans- formation of the mAR coefficients in Eqn. 2.2 into the frequency domain via the Fourier transform A f = I − P ∑ p=1 Ape−i2pi f p, (2.6) where A f is a D×D matrix denoting the Fourier transform of time-domain mAR coefficients Ap, p= 1, . . . ,P, at frequency f with elements A f (i, j), i, j= 1, . . . ,D and I is a D×D identity matrix. The estimate of PDC from the node j to the node i is defined as pii← j( f ) = A f (i, j)√ ∑Dm=1 |A f (m, j)|2 . (2.7) PDC takes on a value between 0 and 1. It measures the relative interaction strengths with respect to a given source signal. The PDC pii← j( f ) describes the linear pairwise related- ness between nodes i and j as a function of frequency after discounting the effect of other simultaneously observed series. To infer the connectivity patterns based on PDC, an important step is to assess the signif- icance level of the estimates. However, as mentioned by Baccala and Sameshima in [13], the theoretical distribution of the PDC measure is generally unknown and, therefore, determina- tion of the estimate significance still remains an open question. Various different techniques have been proposed to address this issue. In this work, we adopt the surrogate data approach described in [19]. First, the signals are transformed by Fourier transform to the frequency do- main, and then a random variable, which is generated uniformly in the range [0,2pi), is added to the original phase of each signal. The surrogate data are obtained by performing inverse Fourier transform on the “phase-randomized” data. By construction, the resulting surrogate data will have the same power spectrum as the original data, but the causal relations between signals are destroyed. This procedure is repeated many times (500 times in this work), and the PDC values derived from these surrogate data form a null distribution. The PDC value from 22 the original signals is compared with the null distribution to obtain a p-value. If the p-value is smaller than a pre-defined significance threshold (e.g. 0.05), a connection is assumed being present. 2.3 Results In this section, we will first examine the PDC results produced by sparse mAR models and by regular mAR models when fitted to simulated data. Next, we will present a case study where the proposed sparse mAR-based PDC is applied to real EEG data collected from an motor task where visual stimuli were presented in a virtual environment. 2.3.1 Simulated Data To compare the PDC estimates produced by sparse mAR models and by regular mAR models, an 18-channel, second-order (i.e. P = 2) mAR system is simulated. The additive noise term e(t) is uncorrelated, white Gaussian noise with zero-mean and unit variance, and the data length is set to be 2,000. Note that, to illustrate the effects of estimation errors in mAR coefficients on PDC values, such a data length is purposely chosen as to be comparable with the length in real EEG data. The mAR coefficient matrices are 18×18 sparse matrices where the values of the non-zero elements are chosen similarly as in [13], and they are given by A1(2,13) = 0.95 √ 2, A2(1,1) =−0.9025, A1(2,1) =−0.5, A2(3,2) = 0.4, A1(10,13) =−0.5, A2(16,18) =−0.2, A1(13,14) = 0.25 √ 2, A2(4,16) = 0.7, A1(10,2) = 0.25 √ 2, A1(5,4) =−0.25 √ 2, A1(5,5) = 0.25 √ 2. First, to compare the performance of regular mAR models and sparse mAR models in de- tecting connectivity between channels, a second-order regular mAR model and a sparse mAR model are separately fitted to the simulated data. Their respective PDC networks including 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1516 17 18 MAR Network (a) PDC Network based on 2nd Order mAR Model 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1516 17 18 Sparse MAR Network (b) PDC Network based on 2nd Order Sparse mAR Model Figure 2.2: PDC networks of simulated data. Only significant connections (p < 0.01) are shown. only significant connections with p < 0.01 (based on 500 surrogate datasets) are shown in Fig. 2.2. The sparse mAR model correctly identifies all causal relations, whereas in the case of the regular mAR model, several connections are falsely identified. Since the model order is typically unknown in practical problems, we next investigate the effects of the model order (i.e., choice of P in Eqn. 2.2) on PDC. Regular and sparse mAR models of order 2, 5, 7 and 10 are fitted to the simulated data. The resulting PDC curves of two selected connections are plotted in Fig. 2.3. We first present an example of no causal relation by looking at the PDC from channel 2 to channel 18 in Fig. 2.3a. Given that 24 A1(18,2) = A2(18,2) = 0, the true PDC is pi18←2( f ) = 0, indicating the absence of direct causal influence from channel 2 to channel 18. When the PDC is calculated based on Ap given by the sparse mAR models, regardless of the choice of model order, the estimated PDC curves (lines with crosses as markers) closely approximate the true PDC (line with squares) at all frequencies, except that when P = 10, the estimated PDC is slightly greater than zero (pi18←2( f ) ∼ 0.01). On the other hand, the estimated PDC given by regular mAR (lines with circles) noticeably deviates from the true PDC and the deviation increases as the model order increases. When the 0.01 significance threshold given by the surrogate data is applied, the PDC estimates given by the regular mAR model, regardless the model order, are being considered statistically significant for most frequencies, even though there is no causal influence from Y2 to Y18 in the true model. Similar observations are also noted in other connections. This ultimately results in falsely identified connectivity patterns. An example of non-zero connectivity that we look at is from channel 13 to channel 2. In this case, as shown in Fig. 2.3b, the true PDC pi2←13 takes on the value of 0.769 at all frequencies. The sparse mAR-based PDC estimates provide very close approximation to the true values and the deviations are around 0.01 for all model orders. In the case of regular mAR, similar to the connection from Y18 to Y12, the PDC estimates differ from the true PDC considerably. One observation worth pointing out is that, when P = 10, the traditional PDC estimate fluctuates substantially across considered frequencies. This artifact can be attributed to the estimation errors incurred when highly complex mAR models are trained by limited data points, and it may lead to misleading interpretations of spectral properties of the under- lying system. 2.3.2 Case Study: EEG Analysis Next, we apply the proposed technique to the EEG data collected from a virtual reality exper- iment. In this experiment, 7 healthy subjects were recruited and asked to respond to visual stimuli in a computer-simulated virtual environment. The stimuli consists of 150 virtual balls that were sequentially launched and loomed directly towards the subject. Among these 150 25 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 Frequency PD C PDC for Y18 <− Y2 True PDC P=2 Regular mAR P=2 Sparse mAR P=5 Regular mAR P=5 Sparse mAR P=7 Regular mAR P=7 Sparse mAR P=10 Regular mAR P=10 Sparse mAR Sig. Threshold of 0.01 (a) PDC for Y18← Y2 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.66 0.68 0.7 0.72 0.74 0.76 0.78 Frequency PD C PDC for Y2 <− Y13 True PDC P=2 Regular mAR P=2 Sparse mAR P=5 Regular mAR P=5 Sparse mAR P=7 Regular mAR P=7 Sparse mAR P=10 Regular mAR P=10 Sparse mAR (b) PDC for Y2← Y13 Figure 2.3: PDC of selected connections. In (a), the line with starts represents the sig- nificance threshold of p = 0.01. Note that most parts of the regular mAR-based PDC curves are above the significance threshold even though there is no causal influence from Y2 to Y18 26 balls, 50% of the balls were distracter balls and the other 50% are target balls. The distracter balls remained white in color during the course of their trajectory whereas the target balls were initially launched as distracter balls (i.e. white in color) and after approximately 1.5 seconds, the target balls revealed itself by changing color from white to blue. Subjects were instructed to ignore distracter balls and reach out to block target balls with a virtual paddle as soon as they turned blue. To record the brain activity, subjects were fitted with an EEG cap with 18 active chan- nels placed according to the international 10-20 placement systems. The placement of the electrodes is shown in Fig. 2.4, and they are, starting from the top left corner, FP1, FP2, F7, F3, FZ, F4, F8, T7, C3, CZ, C4, T8, P7, P3, PZ, P4, O1 and O2. Two additional electrodes were placed above and below the eye for detecting eye-movement artifacts. EEG signals were recorded using the SynAmps2 amplifier and Scan4 software (NeuroScanr, Compumedicsr, Texas, USA). The data were digitally sampled at 1000Hz and later downsampled to 250Hz. To eliminate eye artifacts, independent component analysis (ICA) was applied to the EEG recordings using the procedure described by Jung et al. in [48]. The denoised data were bandpassed between 5-100Hz. For the analysis, only EEG signals collected during the course of successfully blocked target balls are used. Moreover, for each target ball, the EEG signals are segmented into two parts: pre-reveal period, which is between the time the target was launched and the time the target revealed itself, and post-reveal period, which is between the time the target re- vealed itself and the time the target was blocked by the subject. Third order regular mAR models and third order sparse mAR models are separately fitted to the pre-reveal period and post-reveal period. PDC values are calculated based on the estimated mAR coefficients. To study the effects of ball launch, PDC values from pre-reveal period is subtracted from pre- reveal period. Let4PDCregular and4PDCsparse denote the PDC differences given by regular mAR and sparse mAR, respectively. Student t-tests are separately applied to 4PDCregular and 4PDCsparse, and the significant connections (p < 0.05) are shown in Fig. 2.4. The 4PDCsparse results are consistent with demonstrating frontal gamma synchronization in the 27 Significant ∆PDC regular in Gamma Band (a) Regular mAR Model Significant ∆PDC sparse in Gamma Band (b) Sparse mAR Model Figure 2.4: Head plots showing significant connections in Gamma band (36-60 Hz). (a) Significant connections (p < 0.05) as determined by t-test on 4PDCregular (b) Significant connections (p< 0.05) as determined by t-test on4PDCsparse frontal area of cortex in response to complex visual stimuli requiring interpretation [49]. In particular, the4PDCregular and4PDCsparse curves for the connection from CZ to F8, which has the smallest p value in the connectivity network, is shown in Fig. 2.5. The midline of 4PDCsparse curves and the associating standard deviation (shown by the vertical error bars) are consistently below zero, whereas the standard deviation of4PDCregular curves is notice- ably larger and span both positive and negative values, which makes neurological interpreta- tion difficult. 2.4 Conclusion The proposed sparse mAR-based PDC estimates are shown to be more accurate and con- sistent in simulations. In a real EEG study, the results from the sparse mAR-based PDC approach were consistent with prior studies emphasizing frontal gamma oscillations. Our re- sults suggest that when the number of data points are limited as often the case in real word applications, the sparse mAR-based PDC formulation is more suitable for real EEG analysis by accommodating the sparse nature of brain connectivity network and naturally address- ing the concern of model order selection, and that PDC results based on the regular mAR approach should be interpreted with caution. 28 10 20 30 40 50 60 70 80 90 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 Frequency (Hz) ' P D C F8 (ch. 7) <- CZ (ch. 10) F8 CZ Regular mAR Sparse mAR Figure 2.5: Average4PDCregular and4PDCsparse curves across all subjects for the con- nection from CZ to F8. The vertical lines represent the standard deviation. 29 Chapter 3 Altered Directional Connectivity in Parkinson’s Disease Using Sparse PDC 2 3.1 Introduction The Parkinsonian state is associated with abnormal oscillatory activity in the basal gan- glia (BG) and related BG-cortical loops. The local field potential (LFP) recordings in the subthalamic nucleus (STN) and the internal segment of the globus pallidus (GPi) during deep brain stimulation (DBS) surgery for Parkinson’s disease (PD) have revealed the presence of exaggerated synchrony in the 13-30 Hz frequency range [50, 51]. These beta-band os- cillations are suppressed by behaviorally-relevant stimuli [52–54] and levodopa medication [51, 55, 56]. In addition, the degree of levodopa-induced beta-band power suppression in STN LFPs correlates with clinical improvements in both rigidity and bradykinesia [40]. Moreover, direct stimulation of the STN at 20 Hz has been shown to slow motor performance, thus pro- viding further evidence that abnormal BG oscillatory activity in the beta-band may contribute to the slowing of movement typical of PD [57]. In contrast, high frequency (100-180 Hz) DBS of the STN results in suppression of beta-band LFP activity that also correlates with 2A version of this chapter has been published as: G. Tropini, J. Chiang, Z. J. Wang, E. Ty, and M. J. McKeown, “Altered directional connectivity in Parkinson’s disease during performance of a visually guided task,” NeuroImage, vol. 56, no. 4, pp. 2144-2156, Jun. 2011. http://dx.doi.org/10.1016/j.neuroimage.2011.03.018 30 improvements in motor performance [41]. Although prior work on oscillatory activity in PD has emphasized findings from STN recordings, there is increasing recognition that the cortex might be involved in the generation and propagation of abnormal rhythms to a greater degree than has been previously recognized [58]. In fact, abnormal oscillatory activity in PD is not limited to the STN and GPi but rather is a network phenomenon, present throughout the BG-thalamocortical circuitry [51, 56, 59–61]. In addition, high-frequency DBS induces antidromic stimulation of layer V cortical cells in anesthetized healthy and 6-OHDA-lesioned rodents [62] and the size of the antidromic evoked potential correlates with therapeutic improvements from akinesia in awake freely moving rats [63]. Further insight into DBS mechanisms has been obtained via optical deconstruction of the parkinsonian neural circuitry in rats, suggesting that the primary site of DBS action is not the STN, but, in fact, the motor cortex [64]. Furthermore, in humans, single-pulse transcranial magnetic stimulation (TMS) of the primary motor area and supplementary motor area in PD patients at rest and off levodopa medication, has been shown to result in a decrease of beta- band activity in the STN [65]. These results suggest that further exploration into the role of cortical oscillatory activ- ity in PD is warranted. Several previous studies have mainly focused on the investigation of cortical-STN coherence via concurrent scalp EEG (with a limited number of electrodes due to experimental limitations) and depth STN recordings in PD patients [56, 60, 66–69]. However, only a few EEG studies have investigated functional cortico-cortical connectivity in the con- text of abnormal oscillatory activity in PD, again utilizing coherence as a measure of coupling between cortical regions [70, 71]. Other studies have used magnetoencephalogram (MEG) to measure functional connectivity in PD using a synchronization likelihood (SL) technique (a general measure of linear and non-linear temporal correlations between time series) but only investigated the resting state [4]. The present study aims to expand on these previous studies by using a method of analysis that would allow us not only to establish interrelations between EEG nodes or cortical areas, but also to infer causal dependencies between them. PDC methods appear particularly well-suited to investigation of cortical oscillatory activity 31 from the EEG [72–74] as they allow the measurement of the strength of directional connec- tivity between two nodes or regions of interest, (and hence infer a measure of causality) while at the same time discounting the effect of all other neighboring nodes [13]. Since PDC-based methods take into account the influence of other nodes when assessing the connectivity be- tween any given pair of signals, they are also more robust to volume conduction effects that affect standard coherence techniques [75–77]. In order to assess EEG-based cortical connectivity in the context of modulation of ab- normal brain rhythms in PD, we investigated movement preparation and execution during a choice reaction task. This type of paradigm has been commonly used in studies investi- gating LFP activity in the BG, and has been shown to modulate beta-band activity in STN recordings [52–54, 60, 67]. In contrast to prior EEG studies on abnormal oscillatory activity in PD [70, 71] we also assessed modulation of brain rhythms in age-matched normal con- trols, who demonstrated robust reductions in beta-band connectivity between occipital and left sensorimotor regions, likely reflecting the visually-guided nature of the task. 3.2 Methods 3.2.1 Overview Sparse multivariate autoregressive-based PDC (sparse mAR-PDC) connectivity analysis in- cluded a visualization of the temporal dynamics of connectivity during each experimental trial interval, and a quantitative, frequency-dependent analysis of connectivity during both the preparatory and movement phases of the motor task. In addition to connectivity analysis, we also utilized spectral analysis in order to both assess power changes in the regions of in- terest, and to compare our techniques with more traditional methods of analysis. All analyses were conducted using the Matlab Signal Processing Toolbox, and custom Matlab scripts. 32 3.2.2 Sparse mAR-PDC Connectivity Analysis Directional connectivity between pairs of electrode regions (5× 5 connections, minus 5 non- meaningful auto-connections) was determined by computing the PDC connectivity spectrum, based on a 5th order sparse mAR model for each individual trial . Briefly, sparse mAR- PDC involves first fitting of a P-th order mAR model to the data. The sparsity constraint on the mAR model is obtained via penalized regression, which jointly limits the number of parameters to be estimated and reduces computational demands. Once the sparse mAR coefficients are determined, they are transformed into the frequency domain via the Fourier transform, and the estimate of the PDC between two nodes (or regions of nodes) is calculated (refer to Section 2.2 for a detailed description). The PDC takes a value between 0 and 1 and describes the pair-wise strength of directional interaction between two nodes as a function of frequency, while at the same time discounting the effect of other simultaneously active nodes. One major assumption of the mAR model is stationarity of modeled signals. In order to capture the temporal evolution of connectivity patterns during dynamic movements, we used an overlapping time window and the PDC was calculated using the time samples within each window. This technique produced a PDC spectrogram or so-called “PDCogram”, which forms a three dimensional matrix of time, frequency and connectivity strength. The optimal window length and amount of overlapping between successive windows are determined based on the desired temporal and spectral resolutions. 3.2.3 Spectral Analysis For each subject group (Controls, PD-OFF, PD-ON) and each subject, trials were analyzed separately by computing short-time Fourier transform spectrograms (‘spectrogram’ function in Matlab) over the 1-50 Hz frequency range for each individual trial and for the baseline period. As for the sparse mAR-PDC connectivity analysis, trials were truncated to the shortest trial length to ensure consistency. The window size used was 125 samples and the window shift was 16 samples. For each trial and each electrode region, the mean baseline log power spectrum was subtracted from the task spectral estimate to obtain normalized values. Again, 33 for visualization, individual trial spectrograms for each electrode region were averaged within each subject and then averaged across all subjects within one group to obtain group values. 3.3 Experimental Data 3.3.1 Subjects Ten volunteers with clinically diagnosed PD participated in our study (two women, eight men, mean age 60.50 ± 9.71 years, one left-handed). The study was approved by the University of British Columbia Clinical Research Ethics Board (CREB), and all subjects gave written informed consent prior to participation. All patients were recruited from the Movement Dis- orders Clinic at the University of British Columbia, and had mild to moderately severe PD (Hoehn and Yahr Stage 1-3); their clinical details are summarized in Table 3.1. Exclusion criteria included atypical Parkinsonism, presence of other neurological conditions, and visual conditions requiring correction beyond the use of eyeglasses or contact lenses. All patients were on levodopa medication. We also recruited ten healthy aged-matched volunteers (five women, five men, mean age 59.60± 8.96 years, one left-handed) without active neurological disorders. All patients stopped their L-dopa medication overnight for a minimum of 12hrs, and any dopamine agonist medication for at least 18hrs. The mean Unified Parkinson’s Disease Rating Scale (UPDRS) motor score in the off medication (PD-OFF) state was 26.00 ± 8.15. After completing the task in the PD-OFF state, patients were given their usual morning dose of levodopa (180± 58.69 mg). After an interval of approximately 1hr, introduced to ensure that the medication reached peak dose, subjects repeated the task while on medication (PD-ON). Control subjects performed the task only once. However, we deemed practice effects to have reached saturation in both control and PD subjects based on a sufficient number of practice trials, as determined in previous work from our group [78–80]. 34 Table 3.1: Summary of patient details. Subject Disease Age Gender Handedness Motor UPDRS Levodopa Duration Part III (OFF) Daily (yrs) Dose (mg) P004 4 64 M R 30 600 P005 4 44 M R 16 300 P006 5 75 M R 35 400 P007 8 66 M R 25 450 P008 12 62 M R 17 800 P010 9 64 M R 19 800 P011 16 58 F R 22 900 P012 7 55 M R 31 500 P014 4 47 F L 24 400 P016 16 70 M R 41 950 AVG 8.50 60.5 26.00 610.00 STDEV 4.72 9.71 8.15 234.28 3.3.2 Behavioral Paradigm Subjects performed a visually guided choice reaction task using a joystick while sitting 1.5 meters away from a large screen. The task was modified from a paradigm previously shown to modulate beta-band activity in STN LFP recordings [52]. Each trial began with the pre- sentation of a small fixation cross at the centre of the screen. After 250 ms, four circular targets appeared (non-predictive cue) around the central cross (Fig. 3.1). After 2500 ms, a pseudorandomly selected target turned yellow (Go cue). Subjects were instructed to use the joystick to guide a white cursor over the target as fast as possible (within 2000 ms of the Go cue). Once the target was reached, it turned green for 1000 ms, and during this time subjects returned the joystick to its resting position. Lastly, the fixation cross was again presented for 1500 ms. 35 Figure 3.1: Joystick task schematics. “Cue” indicates the non-predictive cue. Each trial was defined as the time from the initial presentation of the central fixation cross, to the end of the 1500 ms interval (last panel on the right). We utilized 68 trials, with an equal number of targets in each direction. Target timing and direction were pseudo-randomized but the same sequence of trials was presented to each subject. Subjects underwent both basic joystick training and specific training on the task (12 trials, as determined in pilot trials). If the subjects moved the joystick prematurely, or failed to reach the target on time, the trial was removed from subsequent analysis. Reaction times (i.e., the start of a movement, as determined by change in initial joystick position) and response times (i.e., the time for subjects to reach a target) were recorded. Rest data were also collected while subjects were fixating on the central cross for 1 minute at the beginning of all the trials. The visual stimuli presented during the task were designed in Matlab (The Mathworks Inc., Massachusetts), using the Psychophysics Toolbox Version 3 [81]. 3.3.3 Data Collection and Analysis Subjects were fitted with an EEG cap (Synamps2 Quik-Cap, Neuroscan, Compumedics, Texas USA) with 19 active channels using the international 10-20 placement system, referenced to the mastoids. Artifacts due to eye movements were recorded by surface electrodes placed 36 above and below the eyes (Xltek, Ontario, Canada). EEG data were recorded at 1000 Hz and stored via the Synamps2 amplifier system and Scan4.3 software (Neuroscan, Compumedics, Texas USA). Alignment of EEG data to task-related events was performed by sending Matlab- coded TTL pulses via the parallel port from the stimulus computer to the Synamps2 system. TTL pulses were then recorded as event-specific timestamps by the Scan 4.3 software. Data were later down-sampled to 250 Hz, and eye and EMG-derived artifacts were re- moved via the Automatic Artifact Removal (AAR) toolbox v1.3 [82] of the EEGLab open source Matlab Toolbox [83]. The denoised data were then bandpassed at 1-100 Hz, normal- ized to unit variance, de-trended (‘detrend’, Matlab) and averaged over five electrode regions. (Fig. 3.2; Fronto-Central (FCentral): Fp1, Fp2, F3, F4, F7, F8, Fz; Left Sensorimotor (LSm): C3, P3, T7, P7; Right Sensorimotor (RSm): C4, P4, T8, P8; Central: Cz, Pz and Occipital: O1, O2). Fp1/Fp2, as well as F3/F4 and F7/F8 have all been shown to correspond to frontal areas, while C3/C4 correspond to Brodmann area 4, and P3/P4 to Brodmann area 7 [84] and were thus grouped together with parietal (P7/P8) and temporal (T7/T8) electrodes as ‘sensorimo- tor’ areas for simplicity and ease of computation. O1 and O2 correspond to Brodmann area 17, the primary visual cortex [84]. Last, Cz and Pz correspond to central and parietal midline areas [85]. For analysis purposes, data were segmented into individual trials. A trial was defined as the time from the initial presentation of the fixation cross to the end of the 1500 ms interval (Fig. 3.1). In order to ensure that trials were all the same length, they were truncated to the shortest trial length (5816 ms). Baseline data were obtained from the rest recordings performed at the beginning of the experimental task, while subjects where fixating on the cross presented at the centre of the screen. For each trial and each electrode region, the mean baseline PDC was subtracted from the task PDC in order to produce baseline-normalized connectivity values. The window size for the PDC connectivity spectrogram was 120 samples (480 ms), the window overlapping was 16 samples (64 ms) and the PDC was computed over the 1-50 Hz frequency range (with 1 37 Fp1 Fp2 Fz Cz Pz F7 F3 F4 F8 T7 P7 C3 P3 C4 T8 P4 P8 O1 O2 Fronto-Central Left Sensorimotor Right Sensorimotor Central Occipital Figure 3.2: Headplot of electrode regions. Electrode placement is based on the 10-20 system. Hz resolution). For visualization purposes, individual trial PDCograms were averaged within each subject and then averaged across all subjects within each group, for each direction of connectivity between all possible pairs of electrode regions. Each PDCogram data point thus provides information on the strength of directional connectivity with respect to rest, at a specific frequency and time, for a pair of nodes or electrode regions. 3.3.4 Movement Preparation (PRE) and Movement (MOV) Phases We quantitatively assessed the results of the temporal dynamics analysis by calculating the average power spectrum and average PDC connectivity spectrum with respect to rest during both the preparatory and movement phases of the motor task. The preparatory phase of movement (‘PRE’) was defined as the portion of the interval preceding the Go cue (0 to 2750 ms), and the movement (‘MOV’) phase was defined as the part of interval immediately following the Go cue and lasting until the average response time for all groups (2750 ms 38 to 3590 ms). Response time was defined as the time the subject had reached the target by using the joystick to guide the white cursor. We utilized response time as opposed to reaction time in order to include the entire extent of the joystick-based movement. The average power spectra and the average PDC connectivity spectra over the 1-50 Hz frequency range were extracted from the temporal dynamics data for PRE and MOV for each individual subject, and mean group values were then computed. Error bars were obtained via leave-one-out cross-validation. To determine the differences in the band-power and connectivity modulations between subject groups, individual trial power and PDC spectral values were averaged over the afore- mentioned frequency bands and one-way analysis of variance (ANOVA) tests (nested design with “subject group” as between factor and “subject” as within factor) were performed be- tween every pair of subject groups (Controls vs. PD-OFF, Controls vs. PD-ON and PD-OFF vs. PD-ON). A separate ANOVA test was performed for each region and connection pair. Significance was established with false discovery rate (FDR)-corrected value of 0.0051 for band-power (initial α = 0.01, corrected for 5 regions, 5 frequency bands and 3 group pairs) , and 0.0050 for PDC (initial α = 0.01, corrected for 20 connections, 5 frequency bands, and 3 group pairs). To assess the clinical significance of the PDC values, individual subject PDC connectivity spectra values for PD-OFF were then averaged over different frequency bands (theta: 4-7 Hz, alpha: 8-12 Hz, beta: 13-30 Hz and low-gamma: 31-50 Hz) and correlated with UPDRS III OFF motor scores, using a robust regression (‘robustfit’, Matlab). Analysis was performed for each connection pair separately. Significance was established with a FDR-corrected alpha level of 0.025 (initial α = 0.05, corrected for 20 connections, 4 frequency bands, and 2 task intervals). 39 3.4 Results 3.4.1 Behavioral Measurements PD-OFF subjects showed increased response times as compared to both PD-ON and age- matched controls (Controls: 0.76 ± 0.15s, PD-OFF: 0.92 ± 0.11s, PD-ON: 0.86 ± 0.13s). Response times for PD-OFF and PD-ON were significantly different from one another when the effect of the medication was examined (Paired t-test, p = 0.035), but only PD-OFF was significantly different from control subjects (Student’s t-test, PD-OFF: p = 0.020 and PD- ON: p = 0.155). A similar trend was found when examining reaction times (Controls: 0.47 ± 0.06s, PD-OFF: 0.51 ± 0.05s, PD-ON: 0.48 ± 0.04s). Again, reaction times for PD-OFF and PD-ON were significantly different from one another (p = 0.014), but not from control subjects (p= 0.206 and p= 0.973, respectively). 3.4.2 Frequency Spectra In normal subjects, spectrograms showed a brief decrease in power relative to baseline in the 10-25 Hz range following the non-predictive cue (Fig. 3.3). Power suppression also occurred approximately 1 s before the Go cue, particularly in the FCentral, LSm and Occipital regions. Following the onset of the imperative cue, a large power suppression in the 13-30Hz range was then observed across all regions. This power suppression lasted until approximately 0.5s following the average response time for each group, presumably due to the fact that sub- jects were instructed to return the joystick to the center position following target selection. Power rebound in the 20-30 Hz range was then observed in all regions approximately 1s after response time, although its magnitude was less pronounced in the Occipital region. Power in- creases at lower frequencies (1-10Hz, and up to 15Hz for RSm, Central and Occipital regions) were observed in all regions following the non-predictive cue. Although the power in these frequency ranges renormalized after about 1s from the beginning of the task, a subsequent re-increase was observed during movement in all regions. In contrast, PD patients OFF medication showed a reduced power suppression in the 40 10-25 Hz range following the non-predictive cue and immediately before the imperative cue (Fig. 3.3). Similarly, power suppression in the 13-30 Hz range was blunted during movement. Lower-band (1-10 Hz) power modulation was also blunted in PD-OFF patients both following the non-predictive cue and during movement, particularly in the Central region. An increased power rebound was observed in the 15-30 Hz range in PD-OFF subjects as compared to control subjects, particularly in the side ipsilateral to movement (RSm) and in the Occipital region. Levodopa medication (PD-ON) partially restored beta-band suppression before and dur- ing movement, albeit not to the levels observed in controls. Power suppression in the 10-25 Hz range following the non-predictive cue was also partially restored. In contrast, PD-ON subjects demonstrated blunted power modulation in the 1-10 Hz both following the non- predictive cue and in correspondence to movement. The medication also reduced the power rebound following movement. The PRE (0-2750 ms) and MOV (2750-3590 ms) intervals did reveal task-dependent differences, while only modest region-dependent differences were observed for all groups (Fig. 3.4). Overall, fewer power differences with respect to rest were observed in the PRE phase as opposed to MOV. In the PRE phase, power was generally reduced with respect to rest with the exception of the right RSm and central regions in normal subjects below 10 Hz. Impairment of beta suppression was seen in PD OFF in LSm, Rsm, and Occipital regions. In the gamma band, reduced power suppression was observed in both PD groups in the FCentral and Occipital regions, while enhanced power suppression in PD ON and OFF was observed in the LSm and Central regions. In the MOV phase, power was increased with respect to rest in Normals < 8 Hz, while a marked decrease in power was observed between 10-35 Hz for all groups. In PD-OFF subjects, beta-band suppression and alpha-band augmentation was significantly less than that compared to controls. While PD-ON subjects had beta-band suppression between that of Normals and PD-OFF subjects, their alpha-band augmentation was in fact worse. In the gamma band, power modulation patterns mirrored those observed in the PRE phase, except 41 Time (s) Mean Reaction Time Mean Response Time F re q u en cy ( H z) Figure 3.3: Spectrograms for all groups. The five electrode regions (FCentral, LSm, RSm, Central and Occipital) are shown on the left. Indicated are power values with respect to rest, in dB. The vertical red bars indicate average response times for each group, while the vertical white bars indicate the average reaction times. “Cue” and “Go” refer to the non-predictive cue and the imperative cue shown in Fig. 3.1. 42 in the LSm region where no significant differences between groups were found. 3.4.3 Connectivity Analysis Overall, sparse mAR-PDC based analysis revealed several asymmetries in the bi-directional patterns of connectivity between the 5 regions of interest, as well as a differential involvement of each region in the temporal dynamics. Connections that appeared to be particularly dis- rupted in PD-OFF subjects were Central↔ FCentral, LSm↔ Occipital (Fig. 3.5), as well as RSm↔ Occipital and Central↔ Occipital (Fig. 3.6). Overall, PDCograms for PD-OFF sub- jects showed dramatic changes in connectivity as compared to controls, while PDCograms for PD-ON subjects revealed that levodopa medication failed to completely restore normal pat- terns of modulation, and that improvement was connection-dependent. In particular, greater medication-dependent improvements were observed for LSm↔ Occipital (Fig. 3.5, bottom panel), and partially for Central→ FCentral (Fig. 3.5, top panel). The PDC-based quantitative analysis of the PRE and MOV intervals aided further identi- fication of several abnormalities in the connectivity modulation by PD subjects. In the PRE phase (Fig. 3.7), PD-OFF subjects showed disrupted connectivity with respect to rest in sev- eral connection pairs. Consistent with that visualized by the PDCograms, information flow was most disrupted in the following connections: Central ↔ FCentral, Occipital → LSm, Occipital↔ RSm, Central→ Occipital. Of these, only the Central→ FCentral connection showed an almost complete return to normal modulation patterns following levodopa admin- istration (PD-ON). Levodopa however further deteriorated connectivity in LSm→ FCentral, LSm→ Central and Occipital→ Central. In order to determine whether eye movements affected the PDC, we extracted the onset of eye movements from the EEG during the artifact-removal processing step. We then assessed the PDC time locked to the eye movements. Specifically, the PDC 320 ms before and 320 ms after the onset of eye movements, when integrated over all conditions, did not show any significant differences (not shown). Correlation of OFF UPDRS Motor Scores with PRE PDC values for each of the PD-OFF 43 0 10 20 30 40 -4 -3 -2 -1 0 1 2 MOV Power Spectra - Occipital0 10 20 30 40 -4 -3 -2 -1 0 1 2 MOV Power Spectra - Central0 10 20 30 40 -4 -3 -2 -1 0 1 2 MOV Power Spectra - R sm 0 10 20 30 40 -4 -3 -2 -1 0 1 2 MOV Power Spectra - L sm 0 10 20 30 40 -4 -3 -2 -1 0 1 2 Frequency (Hz) P o w e r (d B ) MOV Power Spectra - FC 0 10 20 30 40 -4 -3 -2 -1 0 1 2 PRE Power Spectra - Occipital0 10 20 30 40 -4 -3 -2 -1 0 1 2 PRE Power Spectra - Central0 10 20 30 40 -4 -3 -2 -1 0 1 2 PRE Power Spectra - R sm 0 10 20 30 40 -4 -3 -2 -1 0 1 2 PRE Power Spectra - L sm 0 10 20 30 40 -4 -3 -2 -1 0 1 2 Frequency (Hz) P o w e r (d B ) PRE Power Spectra - FC PRE MOV 0 50 FCentral-->L Sm 0 50 FCentral-->R Sm 0 50 FCentral-->Central 0 50 FCentral-->Occipital 0 50 Frequency (Hz) P D C V a lu e s L Sm-->FCentral 0 50 L S -->R Sm 0 50 L Sm-->Central 0 50 L Sm-->Occipital 0 50 R Sm-->FCentral 0 50 R Sm-->L Sm 0 50 R Sm-->Central 0 50 R Sm-->Occipital 0 50 Central-->FCentral 0 50 Central-->L Sm 0 50 Central-->R Sm 0 50 Central-->Occipital 0 50 Occipital-->FCentral 0 50 Occipital-->L Sm 0 50 Occipital-->R Sm 0 50 Occipital-->Central N PDOFF PDON Figure 3.4: Spectra during PRE and MOV phases for the five regions. An ANOVA test was used to assess significant differences in theta, alpha, beta, combined alpha- beta and low-gamma bands. The PRE phase (left column) and MOV phases (left column) are shown. The triangle indicates significant differences between PD OFF and Normal groups, the star indicates significant differences between PD ON and Normal groups, the circle indicates significant differences between PD ON and PD OFF groups (all at significance of p< 0.01). 44 patients (excluding P016, an outlier with UPRS III Score=41) revealed that PDC values in the Central→ FCentral direction of information flow were positively correlated with motor clinical scores in the beta and low-gamma frequency ranges (beta: r = 0.818, p = 0.011; low-gamma: r = 0.55, p= 0.00025 a FDR-corrected α level of 0.025). The correlation with UPDRS UPDRS scores was particularly strong with the low-gamma frequency (Fig. 3.9). In the MOV phase (Fig. 3.8), information flow for PD-OFF was again most disrupted in the Central ↔ FCentral, Occipital → LSm, Occipital ↔ RSm and Central → Occipital. In addition, PD-OFF subjects also showed deteriorated modulation of connectivity for the LSm→ FCentral and LSm→ RSm connections, as compared to the PRE phase. Levodopa medication restored connectivity in the Occipital → LSm direction to levels comparable to those of controls, as well as in the Occipital → RSm (above 10 Hz only) and Central → FCentral directions (above 20 Hz only). As in the PRE phase, levodopa further deteriorated connectivity LSm→ FCentral, LSm→ Central and Occipital→ Central. Correlation of UPDRS OFF Motor Scores with MOV PDC values again revealed a sig- nificant positive correlation between connectivity values for the Central→ FCentral connec- tion. Correlations were again significant in the beta and low-gamma frequency ranges, and approached significance in the alpha range (alpha: r = 0.698, p = 0.028; beta: r = 0.930, p = 0.008; low-gamma: r = 0.9748, p = 0.00012 - Fig. 3.9, FDR-corrected α = 0.025). In addition, since 7 PD subjects were male and only 2 were female, we performed a post-hoc cor- relation of subject gender with PD-OFF PDC-values for each connection and each frequency band. No significant correlations were found in either PRE or MOV, with an FDR-corrected alpha of 0.025. 3.5 Discussion We have shown that cortical activity in PD subjects performing a visually guided task is characterized not only by abnormal changes in spectral power dynamics over time, but also, more dramatically, by altered directional connectivity. Moreover, these changes can be non- invasively and robustly detected with scalp-based EEG recordings. We found multifaceted 45 FCentr Centr LSm Occip Figure 3.5: PDCograms of FCentral ↔ Central (top) and LSm ↔ Occipital for all groups (bottom). The top row of each figure shows left to right direction of con- nectivity and the bottom row shows the reverse direction. Colorbars indicate con- nectivity values with respect to rest (arbitrary units, (a.u.)). The red vertical bar indicates the average response time for each group. ‘N’ indicates Normal controls, ‘PD OFF’ patients off-medication and ‘PD ON’ patients on-medication. 46 RSm Occip Centr Occip Figure 3.6: PDCograms of RSm ↔ Occipital (top) and Central ↔ Occipital (bottom) for all groups. The top row of each figure shows left to right direction of connec- tivity and the bottom row shows the reverse direction. Colorbars indicate connec- tivity values with respect to rest (a.u.). The red vertical bar indicates the average response time for each group. ‘N’ indicates Normal controls, ‘PD OFF’ patients off-medication and ‘PD ON’ patients on-medication. 47 0 10 20 30 40 -0.1 0 0.1 FCentral-->L Sm 0 10 20 30 40 -0.1 0 0.1 FCentral-->R Sm 0 10 20 30 40 -0.1 0 0.1 FCentral-->Central 0 10 20 30 40 -0.1 0 0.1 FCentral-->Occipital 0 10 20 30 40 -0.1 0 0.1 Frequency (Hz) P D C V a lu e s L Sm-->FCentral 0 10 20 30 40 -0.1 0 0.1 L Sm-->R Sm 0 10 20 30 40 -0.1 0 0.1 L Sm-->Central 0 10 20 30 40 -0.1 0 0.1 L Sm-->Occipital 0 10 20 30 40 -0.1 0 0.1 R Sm-->FCentral 0 10 20 30 40 -0.1 0 0.1 R Sm-->L Sm 0 10 20 30 40 -0.1 0 0.1 R Sm-->Central 0 10 20 30 40 -0.1 0 0.1 R Sm-->Occipital 0 10 20 30 40 -0.1 0 0.1 Central-->FCentral 0 10 20 30 40 -0.1 0 0.1 Central-->L Sm 0 10 20 30 40 -0.1 0 0.1 Central-->R Sm 0 10 20 30 40 -0.1 0 0.1 Central-->Occipital 0 10 20 30 40 -0.1 0 0.1 Occipital-->FCentral 0 10 20 30 40 -0.1 0 0.1 Occipital-->L Sm 0 10 20 30 40 -0.1 0 0.1 Occipital-->R Sm 0 10 20 30 40 -0.1 0 0.1 Occipital-->Central PRE 0 50 FCentral-->L Sm 0 50 FCentral-->R Sm 0 50 FCentral-->Central 0 50 FCentral-->Occipital 0 50 Frequency (Hz) P D C V a lu e s L Sm-->FCentral 0 50 L S -->R Sm 0 50 L Sm-->Central 0 50 L Sm-->Occipital 0 50 R Sm-->FCentral 0 50 R Sm-->L Sm 0 50 R Sm-->Central 0 50 R Sm-->Occipital 0 50 Central-->FCentral 0 50 Central-->L Sm 0 50 Central-->R Sm 0 50 Central-->Occipital 0 50 Occipital-->FCentral 0 50 Occipital-->L Sm 0 50 Occipital-->R Sm 0 50 Occipital-->Central N PDOFF PDON Figure 3.7: PDC Spectra for the PRE phase of the task. In the analysis, each column element is the source of connectivity flow, while each row element is the sink of the connectivity flow. The regions of interest are indicated in the diagonal elements. Normals are indicated in blue, PD OFF in red and PD ON in black. The triangle indicates significant differences between PD OFF and Normal groups, the star indicates significant differences between PD ON and Normal groups, the circle indicates significant differences between PD ON and PD OFF groups (all at significance of p< 0.01). 48 0 10 20 30 40 -0.1 0 0.1 FCentral-->L Sm 0 10 20 30 40 -0.1 0 0.1 FCentral-->R Sm 0 10 20 30 40 -0.1 0 0.1 FCentral-->Central 0 10 20 30 40 -0.1 0 0.1 FCentral-->Occipital 0 10 20 30 40 -0.1 0 0.1 Frequency (Hz) P D C V a lu e s L Sm-->FCentral 0 10 20 30 40 -0.1 0 0.1 L Sm-->R Sm 0 10 20 30 40 -0.1 0 0.1 L Sm-->Central 0 10 20 30 40 -0.1 0 0.1 L Sm-->Occipital 0 10 20 30 40 -0.1 0 0.1 R Sm-->FCentral 0 10 20 30 40 -0.1 0 0.1 R Sm-->L Sm 0 10 20 30 40 -0.1 0 0.1 R Sm-->Central 0 10 20 30 40 -0.1 0 0.1 R Sm-->Occipital 0 10 20 30 40 -0.1 0 0.1 Central-->FCentral 0 10 20 30 40 -0.1 0 0.1 Central-->L Sm 0 10 20 30 40 -0.1 0 0.1 Central-->R Sm 0 10 20 30 40 -0.1 0 0.1 Central-->Occipital 0 10 20 30 40 -0.1 0 0.1 Occipital-->FCentral 0 10 20 30 40 -0.1 0 0.1 Occipital-->L Sm 0 10 20 30 40 -0.1 0 0.1 Occipital-->R Sm 0 10 20 30 40 -0.1 0 0.1 Occipital-->Central MOV 0 50 FCentral-->L Sm 0 50 FCentral-->R Sm 0 50 FCentral-->Central 0 50 FCentral-->Occipital 0 50 Frequency (Hz) P D C V a lu e s L Sm-->FCentral 0 50 L S -->R Sm 0 50 L Sm-->Central 0 50 L Sm-->Occipital 0 50 R Sm-->FCentral 0 50 R Sm-->L Sm 0 50 R Sm-->Central 0 50 R Sm-->Occipital 0 50 Central-->FCentral 0 50 Central-->L Sm 0 50 Central-->R Sm 0 50 Central-->Occipital 0 50 Occipital-->FCentral 0 50 Occipital-->L Sm 0 50 Occipital-->R Sm 0 50 Occipital-->Central N PDOFF PDON Figure 3.8: PDC Spectra for the MOV phase of the task. In the analysis, each column element is the source of connectivity flow, while each row element is the sink of the connectivity flow. The regions of interest are indicated in the diagonal elements. Normals are indicated in blue, PD OFF in red and PD ON in black. The triangle indicates significant differences between PD OFF and Normal groups, the star indicates significant differences between PD ON and Normal groups, the circle indicates significant differences between PD ON and PD OFF groups (all at significance of p< 0.01). 49 PRE phase r=0.955, p = 2.5 x 10-4 16 18 20 22 24 26 28 30 32 34 36 -0.15 -0.10 -0.05 0.00 0.05 0.10 ∆ G a m m a P D C UPDRS (motor) MOV phase r=0.9748, p = 1.2 x 10-4 -0.15 -0.10 -0.05 0.00 0.05 0.10 ∆ G a m m a P D C Figure 3.9: Correlation between low-gamma band PDC in the Central→ FCentral con- nection. (top) PRE phase. (bottom) MOV phase. alterations in connectivity that were only partially restored - and often exacerbated - by lev- odopa. Most importantly, PDC-based, task and frequency-specific measures of connectivity were found to correlate with UPDRS motor scores in PD-OFF subjects. A vast majority of seminal studies involving abnormal oscillatory activity in PD have involved patients only, largely because LFP recordings are only available in PD subjects un- dergoing DBS surgery. In contrast, our study aimed to assess task-dependent modulation of cortical activity across regions of interest, and to directly compare PD subjects and normal controls using a paradigm that has been previously shown to modulate STN LFP activity in PD. Simply looking at the power spectrum in the different regions of interest was capable of distinguishing between control and PD subjects (Fig. 3.3). In control subjects, spectrograms showed reduced spectral power relative to baseline levels during movement at 10-30 Hz, in agreement with previous findings investigating EEG event-related (de)synchronization in vi- sually guided motor tasks and choice reaction tasks [86, 87]. In contrast, PD-OFF subjects 50 demonstrated impaired spectral power modulation in the 10-25 Hz range both following the non-predictive cue and during movement preparation, and reduced power suppression in the 13-30 Hz band in correspondence to movement. Attenuated preparatory desynchronization in the alpha and beta bands in PD subjects performing a choice reaction task while ON med- ication has been previously reported [88], and studies describing abnormal activity in the beta band in PD subjects are numerous, both using LFP recordings [51, 55, 59, 60] and MEG [4, 89]. In our study, PD-ON subjects had only partially restored normal modulation of power in the 10-30 Hz range, and blunted modulation in the lower frequencies (< 10 Hz). Other work has shown impaired modulation of low gamma-band activity (31-50 Hz) in PD patients performing an auditory choice reaction task [39]. In particular, patterns of event-related syn- chronization (ERS) and event-related desynchronization (ERD) in response to auditory cues were altered in fronto-parietal EEG electrodes in PD subjects [39]. Similarly, a recent study has shown increased resting-state EEG coherence in the gamma-band (30-45 Hz) in frontal electrode pairs in PD patients’ ON medication [90]. In our study, PD patients’ ON and OFF medication showed a similarly impaired modulation of -band power, particularly in the FCen- tral, Central and Occipital regions. Nevertheless, given that the scalp has a low-pass filtering effect and artifacts from cranial muscle activity occur in gamma-band [91], the gamma-band results should perhaps be interpreted with caution, as the MEG may be more appropriate for further investigations. A key finding of the current study is that although there were overall spectral differences between groups, the differences between regions tended to be relatively modest (Fig. 3.3). In contrast, directional connectivity analysis demonstrated distinctive regional differences within groups, suggesting that altered connectivity may be a more sensitive marker of dis- ruptive oscillations in PD. A previous EEG study demonstrated increased cortico-cortical coherence over the 10-35 Hz range in PD patients at rest [71], but coherence analysis is re- stricted to forcing the connectivity between electrodes/regions to be symmetric. In contrast, PDC-based analysis revealed task and region-dependent changes in connectivity that were clearly asymmetric: for example, the connectivity from the LSm region to the Occipital re- 51 gion was different to that from Occipital to LSm (Fig. 3.4, bottom panel). In PD subjects OFF medication, PDC-based analysis revealed extensive connectivity disruptions that were often not restored, and, in some cases were even exacerbated by levodopa medication (Fig. 3.4, top panel, and Fig. 3.5, bottom panel). We found altered connectivity in PD-OFF subjects in posterior regions, specifically in the Occipital → LSm and Occipital ↔ RSm connections. While it is difficult to associate specific EEG electrodes with particular brain areas, we note that a PDC-based method, by simultaneously taking into account all other regions, will tend to be robust to volume conduc- tion, thus allowing some broad interpretations to be made. The occipital electrodes O1 and O2 may relate to Brodmann area 17, the primary visual cortex [84], while electrodes P3, P4 (in the L and RSm regions) overlie Brodmann area 7, involved in visuomotor coordination [84]. Indeed, PD patients have been known to show impairments in the visual system [92–94] as well as visuo-cognitive deficits [95]. There appears to be a pathophysiologic association between BG dysfunction and occipital glucose metabolism reduction in PD, as indicated by a study which showed an inverse correlation between side-to-side finger tapping asymmetries and occipital glucose hypometabolism in PD patients. More recent work revealed occipital and posterior parietal hypoperfusion in PD patients without dementia that correlated with impaired cortical visual processing [96], while other work has shown a correlation of poor performance in visuospatial and visuoperceptual tests with gray matter atrophy in posterior temporal, parietal and occipital regions [97]. Thus, the impairments of posterior connectivity we observed might in part be related to abnormal visual processing during task performance in PD. Further work will be required to determine if alterations in connectivity correspond to specific psychometric indices in PD. While the spectral differences between PD-OFF and controls were most apparent dur- ing the MOV phase, the connectivity analysis detected abnormalities in both preparation and movement in PD. The altered connectivity seen in PD during the movement preparation phase, in the absence of external predictive cues, is perhaps unsurprising given the difficulty that PD subjects demonstrate when performing internally guided movements [98–101] which 52 may selectively recruit basal ganglia circuits [102, 103]. However, PD-OFF subjects exhib- ited increased connectivity contralaterally (Occipital → LSm), and decreased connectivity ipsilaterally (Occipital → RSm) during both movement preparation and execution, a result that might indicate disease-specific changes in visuo-motor processing. PD subjects increas- ingly recruit externally-guided circuits that rely on somatosensory integration, potentially as a compensatory mechanism [80, 99, 102]. The inability of levodopa to restore these connec- tions implies that visuo-motor changes in PD may still exist in PD subjects on medication. The connectivity analysis was also more sensitive than spectral analysis in revealing the complex effects of levodopa. We note that the inclusion of normal subjects in our study al- lowed us to assess whether or not L-dopa medication “improved” or “worsened” connectivity measures by determining if the PD subjects became more or less like normal subjects after medication. While levodopa worsened (or did not improve) spectral power levels below 10 Hz in all five regions, and particularly during movement, the connectivity analysis revealed region-specific changes induced by the medication. During movement preparation, levodopa fully restored connectivity Central→ FCentral, while during movement execution it best re- stored connectivity Occipital→ LSm. In contrast, during both movement and preparation, the medication had a worsening effect on connections with motor components (LSm→ FCentral and LSm→ FCentral). Further work is required to determine if side effects of levodopa ther- apy, which include motor fluctuations and dyskinesias [104, 105], are partially mediated by altered connectivity. Our results indicate that altered connectivity in PD subjects OFF medication correlates with motor UPDRS scores in the beta and alpha-beta frequency ranges during both move- ment preparation and execution, as well as in the alpha frequency range during movement only. Altered modulation in the Central→ FCentral connection might again underlie deficits in visuo-motor processing along a parietal to centro-frontal stream, given that ‘Central’ elec- trodes (Cz, Pz) are located over the central midline in parietal (Pz) and central areas (Cz is located in front of the central sulcus [85]). However, more work is needed to elucidate the role of this connection via specific visuomotor tests, as well as purely motor or visual paradigms, 53 given that our task included both motor and visual components. The finding that PDC values for this connection consistently correlated with clinical measures during both phases of the task also raises the possibility that altered connectivity Central→ FCentral might be relatively robust to a specific task, and may be possibly utilized in the future as a potential biomarker for PD deficits. We chose to average the EEG waveforms from the electrodes in each region of interest, but there may be more optimal methods for identifying the contribution of each electrode to the regions of interest. We explored an ICA-based method to determine optimal weighting of individual electrodes (data not shown), but the results obtained were not significantly different from those obtained here by the simpler averaging method. We note that by averaging the potentials of neighboring channels within a brain region, the results are less sensitive to slight variations in EEG cap placements across subjects, likely reducing intersubject variability. We could have chosen to group the electrodes into more or less than the five groups currently used. Since the PDC must take into account all other connections, we have found that increasing the number of electrode groups (and hence simultaneous connections) to much larger numbers (e.g., 20) results in numerical instability. This is likely due to the fact that we would need a large window size, which may violate the stationarity assumption required for AR model estimation. Already we are looking at relatively short time periods (e.g., the response time), so we do not have room for long window sizes. Furthermore, the averaging technique effectively reduced the number of signals from 20 to 5 and the number of potential connections from 380 to 20, which makes the results much more easily interpretable. For a full 5th order model, using 20 signals would result in 1900 parameters to be estimated. We did use a sparse multivariate autoregressive model to calculate the PDC both to reduce the number of parameters to be estimated, and this in fact may be a more appropriate for EEG analysis as a full-connectivity assumption (i.e., the assumption that each EEG electrode contains activity that is functionally connected to all other nodes). The full-connectivity assumption is at odds with studies demonstrating that the EEG has “small world” network properties. That is, most nodes are not directly connected to one another, but cluster of nodes may communicate via 54 long-range connections [106, 107]. We therefore chose five regions as a reasonable tradeoff between stability of the PDC calculations and ease of interpretability of the results. We focused on the period after the cue to assess the movement phase PDC values. This inevitably combined the time before and after reaction time and may be a potential confound if PDC values were different between these phases and reaction times were different across groups. Ideally, it would have been desirable to divide the response time into the pre-reaction and post-reaction times to calculate the PDC. However, technical factors prevent us from doing so: we are dealing with a 5th order (albeit sparse) mAR model with many parameters to estimate, and we believe we would not get reliable models dividing the overall response time. Furthermore although reaction times were different between PD-OFF and PD-ON, there were no significant differences between PD-OFF and controls and PD-ON and controls. Lastly, due to recruitment limitations, one PD subject was left-handed. We chose to have subjects perform the task with their dominant hand, as the joystick task requires significant dexterity, and would have thus created difficulties in movement execution if subjects were to use their non-dominant hand instead. While subject groups were still directly comparable, as one control subject was also left-handed, it would have been desirable to have greater hand- edness consistency in our subjects. However, the use of leave-one-out validation techniques did not reveal significant changes in our results when excluding the left-handed subjects from the analysis. Our results suggest that the use of sparse mAR-PDC-based analysis might be ideally suited to detect connectivity modulation between cortical regions as a function of frequency. While EEG-based methods do not provide the level of spatial resolution offered by fMRI tech- niques, they allow for excellent temporal resolution and provide information on oscillatory activity in frequency bands of interest. Thus, the use of an EEG analysis method that also pro- vides directional connectivity information -as offered by some fMRI-based techniques [108]- can be beneficial when studying neurological conditions such as PD, that are characterized by abnormal oscillatory activity in cortical and subcortical systems. In the case of PD, infor- mation on the modulation of cortical connectivity is particularly valuable given recent results 55 pointing to a more important role of the cortex in the modulation of BG oscillatory activity than previously recognized [58, 62–64]. Thus, our results suggest that the role of the EEG to monitor abnormal oscillatory ac- tivity in PD might need to be expanded. In fact, non-invasive findings from the EEG could potentially benefit existing therapies for PD or even open new therapeutic avenues. The ef- fects of DBS on cortical connectivity could be for example assessed via EEG measures, and stimulator settings could be subsequently optimized to reduce altered connectivity. Similarly, given the preliminary result that non-invasive TMS can reduce abnormal rhythms in the BG [65], results from the EEG could be potentially used to direct TMS stimulation and optimize restoration of connectivity. Our findings might also increase our understanding of the pathophysiology of PD, par- ticularly in the context of the oscillatory model of the basal ganglia. Frequency-dependent abnormal modulation of connections associated with visual and visuo-motor processing sug- gests that altered processing in these pathways in PD could be due not only to structural ab- normalities or changes in the firing rate of the fibers involved, but potentially also to altered synchronization and desynchronization patterns, although further work is needed to address these speculations. Similarly, the finding that levodopa results in deteriorated modulation of connections that might be associated with motor processing raises the possibility that altered connectivity might be related to the motor side effects of levodopa such as levodopa-induced dyskinesias (LID). 3.6 Conclusions We have shown that PD subjects demonstrate impairments not only in the intrinsic modula- tion of activity in each region of interest -as shown by spectral analysis- but also profound and complex alterations in modulation of directional connectivity that might underlie deficits in visuo-motor processing. In addition, our work shows that levodopa medication does not fully restore connectivity, and in some cases further exacerbates the deficits shown by unmedi- cated subjects. Moreover, changes in connectivity that might underlie visuo-motor processing 56 correlate with motor UPDRS scores in unmedicated patients in frequency bands previously associated with abnormal oscillatory activity in PD (alpha and beta). Our results suggest that the EEG, coupled with appropriate analyses, demonstrates abnor- malities in connectivity in PD subjects performing paradigms previously shown to modulate LFP beta band oscillations. Given the widespread availability of the EEG, the role of the EEG to monitor PD may need to be further expanded. 57 Chapter 4 Prediction of PD Severity Using Sparse PDC Connectivity Strengths 4.1 Introduction Previous works have shown that altered beta-band oscillatory activity exists in PD, and may relate to the clinical feature of bradykinesia. More recently, studies involving functional con- nectivity at rest in PD have found correlations between connectivity and clinical measures. Using an eyes-open paradigm, Silberstein et al. [71] found that increases in EEG-EEG coher- ence correlated positively with disease severity, while reductions in coherence correlated with improvements in motor function [71]. An MEG study employing an eyes-closed paradigm found that severity of parkinsonism was positively associated with SL values (a general mea- sure of linear and non-linear temporal correlations between time series) in the theta and beta bands, while disease duration correlated positively with SL in the high alpha and beta bands [4]. Previous work from our group [109] has also shown that directional PDC connectivity along a parietal to fronto-central stream correlated positively with UPDRS III motor scores in the beta and gamma bands in PD off-medication patients during a visually guided motor task. In contrast, two MEG studies that directly assessed power changes in either both eyes- open and eyes-closed conditions [110] or in an eyes-closed condition only [111] did not find 58 any significant correlations between power modulation and clinical measures. These findings suggest that functional connectivity might be a better predictor of disease-related measures than changes in oscillatory activity alone at individual brain regions. In recent literature, EEG and MEG functional connectivity measures have been proposed as potential biomarkers for a variety of neurological and psychiatric conditions, such as epilepsy, bipolar disorder and depression, and have even shown the potential to be used as markers for genetic differences in brain network organization in healthy individuals at rest. In this study, our goal to investigate the relationship of EEG-derived cortico-cortical con- nectivity measures during movement preparation and execution with clinical measures of motor function in PD subjects. Specifically, we used a sparse mAR-based PDC method to ascertain directional patterns of connectivity between brain regions. A novel feature selec- tion technique is developed to select connection and frequency band combinations that best differentiate PD subjects from normal subjects and predict disease severity of PD. In clinical settings, the UPDRS has been the “gold standard” assessment tool for char- acterizing impairments in PD patients but it is time-consuming and can only be performed by trained medical personnel. Previous studies have demonstrated the possibility of using features from speech recordings to predict PD severity [43, 44]. Our approach of using brain connectivity as predictive features of disease severity can have useful clinical implications in understanding the effects of disease progression in modulating brain connectivity. 4.2 Methods 4.2.1 Feature Selection To identify brain connections that are able to predict PD severity, we model the relationship between PDC values and PD severity as a linear regression: Y = Xβ + e, where Y is the response vector containing PD severity index of all PD-OFF and control subjects (here, we assume that PD severity index is 0 for all control subjects) and X is the predictor matrix with rows corresponding to the subjects and columns containing the PDC values of all possible 59 pairs of brain regions and all frequencies. We further assume a priori that the weighting vector β is sparse, i.e., only a small fraction of elements of β are non-zeroes. To learn a noisy, under-determined linear system which is presumed to be sparse, we incorporated the LASSO (least absolute shrinkage and selection operator) regularizer [112] into the least- square estimate of the regression model as follows: min‖Y −Xβ‖2,subject to ‖β‖1 ≤ λ , (4.1) where λ is the shrinkage running parameter that determines the level of sparsity the solution will be. To solve Eqn. 4.1, we used the least angle regression (LARS) algorithm [113] as it is efficient in finding the LASSO solution from the null model and thus is particularly suitable for our highly under-determined learning problem. To determine the optimal value of λ , we train the regression model on a set of λ values and the goodness-of-fit of each learned model is measured using the coefficient of determination R2 which is defined as R2 = 1− SSerr SStot , (4.2) where SSerr = ∑i(yi− ŷi)2 and SStot = ∑i(yi− ȳ)2. The terms yi and ŷi denote respectively the true and predicted PD severity of subject i and ȳ denotes the mean of true PD severity of all subjects. Assuming the models are sorted in an ascending order based on their associated λ values, the first model that has R2 greater than 0.9 is selected as the optimal model. The non-zero coefficients in this optimal model are further tested against the zero null hypothesis by applying a standard t-test with confidence level at 0.95. 4.2.2 Cross Validation In the conventional approach to predictor development, the data set is separated into a training set for predictor construction and a test set(s) for assessment of the predictive performance. Using the training set, a subset of features is selected from which a regression-based predic- tive model is constructed and applied to the test set. In this study, we assess the predictive 60 ability of the selected features using leave-one-subject-out cross validation procedure which involves leaving out one subject as the test set and using the remaining subjects as the training set. This is repeated such that each subject is used once as the test set. However, in our pre- liminary analysis, we found that the features selected tend to vary with the subjects included in the training set and yield unsatisfactory predictive performance. To circumvent this issue, we propose using a jackknife procedure in which we sequentially leave out one subject at a time from the training set and apply the LASSO-based feature selection technique to the remaining subjects in the training set. This step yields a list of predictive features from each jackknife iteration. The predictive features from all jackknife iterations are pooled together and the frequency of occurrences for each of the selected features is determined. The features that have more than m occurrences (m= 2 in this study) are retained. Based on these features, a linear regression model is constructed (assuming that the i-th subject is left out as the test set): Y traini = X̃ train i βi, (4.3) where the response vector, Y traini , contains the UPDRS values of the subjects in the training set and the predictor matrix, X̃ train i , contains the selected feature variables (with frequency of occurrences≥ 2) of the subjects in the training set. The regression coefficients β̂i is estimated using ‘regress’ function in Matlab. The predicted UPDRS value of the subject in the test set is given by ŷi = x̃iβ̂i, where x̃i is a row vector consisting of the selected feature variables of the test subject. Note that the jackknife-based feature selection procedure and predictor development occurred completely independent of the test set (i.e., only data from the training set is used) and thus cross-validation is unbiased. An illustration of the leave-one-subject-out cross validation procedure is shown in Fig. 4.1. 4.3 Experimental Data Subjects recruited for this study are the same as those recruited for Chapter 3. In this study, subjects were asked to perform a similar visually guided choice reaction task. The task was 61 Subject i Subject 1 Subject 2 Subject i-1 Subject i+1 Subject N Subject 3 Test Set Training Set Subject i Leave-1-out Cross Validation 1 2 N 3 4 … . … . … . 1 2 N 3 4 … . 1 2 N 3 4 … . Jackknife Procedure LASSO LASSO LASSO …. 0 5 10 Connection # O c c u rr e n c e s Selected feature subset i train i train i XY ~ Select corresponding columns of training data to form train iX ~ Regression Analysis îiii xy ̂ ~ˆ Predict UPDRS of Subject i Jackknife Samples Figure 4.1: Illustration of cross validation procedure. In each cross-validation iteration, one subject is left out as the test subjects while the remaining subjects serve as training set. Next, we sequentially leave out one subject (shaded in grey) at a time from the training set and apply LASSO feature selection technique to the remaining subjects (shaded in yellow). The predictive features selected from each jackknife iteration are pooled together to form a histogram. The features that are have more than m occurrences (above the dashed line) are retained. A predictor based on linear regression model is constructed from the training set and used to predict the UPDRS of the subject in the test set. 62 17 Figures Figure 1. Behavioral Paradigm Figure 4.2: Joystick task schematics with predictive warning cue. modified from the one in Chapter 3 to include predictive warning cues as in recent work by Williams et al. [54]. Each trial began with the presentation of a small fixation cross at the centre of the screen. After a variable delay of 2500-3500 ms, introduced to reduce prediction of the timing of the warning cue, four circular targets appeared around the central cross (Fig. 4.2). A predictive cue (yellow arrow) was presented simultaneously for 500 ms. After 2000 ms, a randomly selected target turned yellow (Go cue). Subjects were instructed to use the joystick to guide a white cursor over the target as fast as possible (within 2000 ms of the Go cue). Once the target was reached, it turned green for 1000 ms, and during this time subjects returned the joystick to its resting position. Lastly, the fixation cross was again presented for 1500 ms. Data acquisition and preprocessing followed the same steps as in Section 3.3.3. For analysis purposes, data were segmented into individual trials. A trial was defined as the time from the initial presentation of the fixation cross to the end of the 1500 ms interval (Fig. 4.2). In order to ensure that trials were all the same length, they were truncated to 63 the shortest trial length (5816 ms). Baseline data of each subject were obtained from the rest recordings performed at the beginning of the experimental task, while the subject was fixating on the cross presented at the centre of the screen. Each trial was segmented into the preparatory and movement phases of the motor task. The preparatory phase (‘MOV-PREP’) was defined as the portion of the trial immediately following the end of the presentation of the warning cue to the start of the GO cue (lasting for 2000 ms). The movement phase (‘MOV-EXEC’) was defined as the portion of the trial immediately following the start of the GO cue and lasting until the average response time for all groups (2750 ms to 3590 ms). Response time was defined as the time the subject had reached the target by using the joystick to guide the white cursor. We utilized response time as opposed to reaction in order to include the entire extent of the joystick-based movement. The PDC connectivity between pairs of electrode regions (5 × 5 connections, minus 5 non-meaningful auto-connections) was computed over the 1-50 Hz frequency range based on a 5th order sparse mAR model (refer to Section 2.2 for a detailed description of sparse PDC calculation). The PDC of the MOV-PREP and MOV-EXEC phases are computed separately for each trial. The mean baseline PDC was subtracted from the task PDC in order to produce baseline-normalized connectivity values. The average PDC of each subject was obtained by averaging PDC over trials for each connection. For comparison, the power spectra over the 1-50 Hz frequency range were also computed for MOV-PREP and MOV-EXEC phases of each trial. The average power spectra of each subject were obtained by averaging power spectra over trials for each brain region. The power spectra were then averaged over different frequency bands (theta: 4-7 Hz, alpha:8-12 Hz, beta: 13-30 Hz and low-gamma: 31-50 Hz). 4.4 Results Particular combinations of the cortico-cortical directional connectivity strengths measured by PDC are able to accurately predict UPDRS scores and differentiate PD-OFF subjects from control subjects. The predicted UPDRS scores for the MOV-PREP and MOV-EXEC phases 64 Figure 2 (b) (a) 0 10 20 30 40 -20 -10 0 10 20 30 40 50 UPDRS score P re d ic te d U P D R S u s in g B ra in C o n n e c ti v it y MOV-PREP r = 0.74 (p = 2.6e-004) 0 10 20 30 40 -60 -40 -20 0 20 40 60 UPDRS score P re d ic te d U P D R S u s in g B a n d P o w e r MOV-PREP r = -0.14 (p = 0.5660) Figure 4.3: Prediction performance for the movement preparation phase. (a) Predicted UPDRS based on PDC of selected connections. Note the clear separation between normal and PD subjects. (b) Predicted UPDRS based on average band power. are shown in Fig. 4.3(a) and Fig. 4.4(a), respectively. For the MOV-PRE phase, the features selected by the LASSO approach concentrate in the alpha band (8-12 Hz). The correlation between the true UPDRS scores and the predicted values is 0.74 (p< 0.0001). The predicted UPDRS values of PD-OFF and control subjects are significantly different (p< 0.0001 based on t-test) and yield a misclassification rate of 5.3%. The features selected by LASSO for predicting UPDRS are mainly connections coming out of Central region (Central→ FCentral at 8 Hz, Central→ LSm at 11 Hz, Central→ RSm at 8 and 11 Hz) and connections coming out of RSm area (Rsm→ Central at 8 Hz and RSm→ Occipital at 8 Hz) (see Fig. 4.5). For the MOV-EXEC phase, the features selected by the LASSO approach concentrate in the combined alpha-beta band (8-35 Hz). The correlation between the true UPDRS scores and the predicted value is 0.87 (p< 0.0001). The predicted UPDRS values of PD-OFF and control subjects are significantly different (p = 0.0002 based on t-test) and yield a misclassification rate of 10.5%. Similar connections as those used in MOV-PREP phase were selected for predicting UPDRS for the MOV-EXEC phase, yet wider frequency ranges, encompassing both alpha and beta bands, were selected: Central→ FCentral at 25-28 Hz, Central→ LSm 65 Figure 4 (b) (a) 0 10 20 30 40 -10 0 10 20 30 40 UPDRS score P re d ic te d U P D R S u s in g B ra in C o n n e c ti v it y MOV-EXEC r = 0.87 (p = 1.5e-006) 0 10 20 30 40 -20 -10 0 10 20 30 40 50 UPDRS score P re d ic te d U P D R S u s in g B a n d P o w e r MOV-EXEC r = 0.17 (p = 0.4843) Figure 4.4: Prediction performance for the movement execution phase. (a) Predicted UPDRS based on PDC of selected connections. Note the clear separation between normal and PD subjects. (b) Predicted UPDRS based on average band power. at 15-24 Hz, Central → RSm at 8 and 18 Hz, FCentral → Central at 29-34Hz and RSm → Occipital at 8-10 Hz (see Fig. 4.6). We further performed prediction of the UPDRS scores of PD-ON subjects using connec- tivity from the MOV-EXEC phase, with PD-OFF and control subjects as the training set and PD-ON subjects as the test set. Fig. 4.6 shows the predicted UPDRS scores of subjects in each subject group (control, PD-OFF and PD-ON) where the arrows represent the effects of medication for each PD subject. We see that for most of the PD subjects the predicted UP- DRS values were reduced to a similar level as controls after medication, demonstrating the normalizing effects of levodopa. For comparison, we performed prediction of the UPDRS scores using the average power of all frequency bands and all brain regions (5 bands x 5 regions). The predicted UPDRS scores for the MOV-PREP and MOV-EXEC phases are shown in Fig. 4.3(b) and Fig. 4.4(b), respectively. The average band power is unable to accurately predict the UPDRS scores and the predicted UPDRS values of PD-OFF and control subjects are not significantly different (p= 0.5660 for MOV-PREP and p= 0.4843 for MOV-EXEC) 66 20 20 Figure 4. MOV-PREP phase: features selected by LASSO and their respective number of times being selected in the jackknife procedure. Figure 4.5: Features selected by LASSO for the MOV-PREP phase and their respective number of times being selected in the jackknife procedure. 4.5 Conclusion Previous work has suggested that altered beta-band oscillatory activity exists in PD, and may relate to the clinical feature of bradykinesia. In this work, we proposed a novel feature selec- tion technique and showed that particular combinations of connections/frequency bands were able to clearly differentiate PD subjects from controls and accurately predict clinical UPDRS 67 22 22 Figure 6. MOV-EXEC phase: features selected by LASSO and their respective number of times being selected in the jackknife procedure. Figure 4.6: Features selected by LASSO for the MOV-EXEC phase and their respective number of times being selected in the jackknife procedure. 68 23 23 Figure 7. Predicted UPDRS values of each subject group using connectivity from the MOV-PREP phase. The arrows represent the effect of medication for each PD subject. For most PD subjects, the medication reduced the predicted UPDRS values. Control PD-OFF PD-ON -20 -10 0 10 20 30 40 P re d ic te d U P D R S u s in g B ra in C o n n e c ti v it y MOV-PREP Figure 4.7: Predicted UPDRS values of each subject group using connectivity from the MOV-PREP phase. The arrows represent the effect of medication for each PD subject. For most PD subjects, the medication reduced the predicted UPDRS values. scores in PD subjects. The frequencies selected by the proposed technique were not restricted to the beta band but included both alpha and beta bands, complementing the observation made based on STN recordings. The results suggest that EEG-derived connectivity may be used as an objective measure for evaluating PD deficits. 69 Chapter 5 Sparse Group mAR Model for Assessing Effects of Galvanic Vestibular Stimulation on Brain Connectivity 5.1 Introduction The vestibular system may be considered as a sixth sense [114]; yet, unlike other sensory systems, thalamic and cortical processing is complex, multimodal and widespread. While the parieto-insular vestibular cortex has been described as the “core” vestibular region in non-human primates [115], recent views deviate towards the notion of a highly distributed vestibular network comprising: the temporo-parietal junction, posterior insula, somatosen- sory cortex, posterior parietal cortex, anterior insula, lateral and medial frontal cortices [116]. Due to the widespread thalamocortical projections, our understanding of the spectrum of vestibular influences on cortical and sub-cortical structures remains incomplete. The galvanic vestibular stimulation (GVS) offers significant advantages to investigate the effect of vestibular input on brain function. Without exciting other sensory systems, transcu- taneous delivery of galvanic current to the mastoid processes alters the firing rates of primary vestibular afferents [117]. Direct and precisely controlled activation of the vestibular system 70 avoids unwanted side effects of discomfort and imbalance [118] and has highly facilitated the modern study of balance and dynamic movements [119]. More recently, Yamamoto et al. [120] delivered subthreshold, zero-mean, noisy (i.e., with random fluctuations) GVS in the context of motor tasks to patients with central neurodegenerative disorders, including PD. Patients improved in their motor responsiveness during periods of stimulation, a finding that has been subsequently reproduced [121], although not yet fully explained. The reported motor benefit of noisy GVS in PD is particularly intriguing given that the Parkinsonian state is characterized by abnormal brain oscillations. Neural interactions be- tween the subthalamic nucleus and the external globus pallidus act as one of the primary pacemakers of oscillatory activity around a frequency of 15-30 Hz [61, 122], which can be exaggerated in PD [38] and may serve as the basis of bradykinesia [40]. Since a number of other diseases may also be characterized by abnormal oscillations, including neurogenic pain, tinnitus, and depression [123], the effect of GVS on ongoing – normal and pathological – brain oscillations is an unaddressed issue of further interest. Here, we questioned whether imperceptible, noisy GVS is capable of modulating alpha (8-12 Hz) and beta (12-30 Hz) rhythms. Previously, large EEG artifacts have disrupted the on- going measurement of microvolt-level brain oscillations, a complication now circumvented by: 1) improved EEG amplifier design with high common-mode-rejection ratio and 2) en- hancement of (computationally intensive) artifact rejection analytical techniques such as ICA [124]. Using subsequent spectral density analysis and PDC methods to assess cortico-cortico connectivity [109], we demonstrate that subthreshold, noisy GVS modulates the connectivity and amplitude of ongoing brain rhythms in healthy subjects. 5.2 Methods 5.2.1 Spectral Analysis To assess the relative changes in amplitude within alpha and beta bands, we first computed the power spectral density in each region over each stimulation window (using the “psd” 71 function in Matlab, with nFFT = 1024, window = 512 points, overlap = 256 points). At each stimulation level and each region, we then computed the relative alpha and beta powers by summing the values between 8-12 Hz and 12-20 Hz, respectively, and dividing by the total power over all frequencies. We focused on the alpha and beta values as these demonstrated observed changes with respect to stimulation current. Robust regression statistical analysis (using the “robustfit” function in Matlab) was utilized to determine the correlation between relative alpha and beta powers and stimulation current [125], and significance was established at p values less than 0.05. Specifically, the relative power in each region of interest (ROI) was modeled as: Prelative = β0+β1I, (5.1) where Prelative was the relative power (alpha, beta) at each ROI, β0 was the mean relative power, and β1 was the linear relation between relative power and GVS current intensity (I). We examined the directional connectivity between pairs of electrode ROIs, resulting in 20 comparisons (5 × 5 connections, minus 5 non-meaningful auto-connections), by computing the PDC connectivity spectrum. PDC methods measure the strength of directional connec- tivity between two ROIs while simultaneously discounting the effect of all other nodes [13], and are therefore appropriate for assessing EEG cortical oscillatory activity [73] and tend to be robust to volume conduction effects compared to say, standard coherence techniques [76]. The PDC spectra we used were based on a 5th order sparse mAR model for each individual trial using a Group LASSO (gLASSO) method [126]. The mAR model captures the causal influences between electrode ROIs by representing the current EEG measurements as a lin- ear combination of previous measurements whereas incorporation of gLASSO into the mAR model estimation promotes sparsity in the connectivity pattern at the group level. Specifi- cally, gLASSO enforces the mAR models of the subjects within a group to be structurally identical across subjects, but the connection coefficients are allowed to vary on a subject-by- subject basis. This is done by formulating the group mAR estimation as a constrained least square minimization problem with the penalty term being the sum of `2 norms of the subject- 72 specific coefficients of each connection. Estimates of the mAR coefficients are determined by gLASSO using the CVX optimization package [127], and are consequently transformed into the frequency domain via the Fourier transform. The estimate of the PDC between two electrode ROIs is calculated and normalized, taking values between 0 and 1. 5.2.2 Sparse Group mAR Model Suppose there are D ROIs and the EEG measurements are taken at time 1, . . . , N.We further assume that there are S subjects. Let xsi (t) denote the EEG measurement of subject s of i-th ROI at time t, for i ∈ {1, . . . , D} and s ∈ {1, . . . , S}. Before introducing the group mAR model, we first define the mAR model at the subject-level: xs(t) = P ∑ p=1 Aspx s(t− p)+es(t), (5.2) where xs(t) = [ xs1(t) x s 2(t) . . . x s D(t) ] T is a D-dimensional column vector denoting the EEG measurements of subject s across all D ROIs at time t and P is the mAR model order. The mAR coefficient matrix Asp, p = 1, . . . ,P, is a D×D matrix with elements Asp(i, j) de- scribing the influence of jth ROI on ith ROI after p time points for subject s. The vector es(t) is the innovation noise which is assumed to follow a multivariate white Gaussian process. Putting all subjects together and rewriting in a matrix form, we obtain a linear regression representation of the group mAR model for each target ROI i: x1i, P+1:N x2i, P+1:N ... xSi, P+1:N ︸ ︷︷ ︸ ,yi = X 1 X 2 . . . X S b1i b2i ... bSi + e1i e2i ... eSi ︸ ︷︷ ︸ ,Xbi+ei (5.3) 73 where we define xsi, P+1:N = [ xsi (P+1) x s i (P+2) . . . x s i (N) ]T ∈R(N−P)×1, (5.4) X sp = [ xs1, P+1−p:N−p x s 2, P+1−p:N−p . . . x s D, P+1−p:N−p ] ∈R(N−P)×D, (5.5) X s = [ X s1 X s 2 . . . X s P ] ∈R(N−P)×DP, (5.6) bsi = [ As1(i,1) . . . A s 1(i,D) A s 2(i,1) . . . A s 2(i,D) . . . A s P(i,1) . . . A s P(i,D) ]T ∈RDP×1. (5.7) The gLASSO solution is obtained by formulating the group mAR estimation as a constrained least square minimization problem with the penalty term being the sum of `2 norms of the subject-specific coefficients of each connection: b̂i = argmin‖yi−Xbi‖22+λ D ∑ j=1 P ∑ p=1 ∥∥∥bgLASSOi, j (p)∥∥∥2 (5.8) where λ is the shrinkage parameter which determines how sparse the solution is. The vector bgLASSOi, j (p) = [ A 1 p(i, j) A 2 p(i, j) . . . A S p(i, j) ] T contains the subset of bi corresponding to the connection from of jth ROI to ith ROI at a lag of p time points of all S subjects. The design of the penalty term in Eqn. (5.8) ensures that the coefficients in each sub-vector, bgLASSOi, j (p),∀ i, j ∈ {1, . . . , D} , p ∈ {1, . . . , P}, are either all zeros or all nonzeros. The estimates of the full mAR coefficient matrices, Asp, are obtained by solving Eqn. (5.8) for all i= 1, . . . , D. 74 5.3 Experimental Data 5.3.1 Subjects Ten healthy subjects (five men, five females; mean age ± SD, 37.2 ± 17.7 years) without any reported vestibular, auditory or neurological disorders participated in the primary study. The University of British Columbia Ethics Board approved the study, and all subjects gave written, informed consent prior to participating. 5.3.2 Stimulus GVS was delivered to subjects through carbon rubber electrodes (17 cm2) in a binaural, bipo- lar fashion [117]. Electrodes were placed over the mastoid processes and were coated with electrode gel (TAC-Gel, Pharmaceutical Innovations Inc., NJ, United States) to optimize con- ductivity and adhesiveness. Digital signals were first generated on a PC computer with Lab- view software (National Instruments, TX, United States) and then converted to analog sig- nals via a NI USB-6221 BNC digital acquisition module (DAQ; National Instruments). The analog command voltage signals were subsequently passed to a constant current stimulator (Digitimer DS5 Isolated Bipolar Constant Current Stimulator, Digitimer Ltd., Hertfordshire, United Kingdom), which was connected to the stimulating electrodes. The delivered stimu- latory waveform was digitized by the DAQ in conjunction with Labview software (National Instruments) for additional analysis. Stimulation signals were zero-mean, linearly detrended, noisy currents with a frequency range of 0.1-10 Hz and a 1/f-type power spectrum (pink noise), sampled at 1 kHz (Fig. 5.1). 5.3.3 Primary Study Protocol After initial set up and obtaining proper consent, subjects underwent a threshold test to deter- mine the individual sensory threshold to vestibular stimulation. Starting from a basal current level, GVS was delivered to subjects for 20-second periods with step-wise increases in cur- rent intensity. The sensory threshold level was defined as being reached when the subject 75 Figure 5.1: Characteristics of stimuli delivered. A representative output signal delivered to a subject through a constant current stimulator. Top: Frequency distribution of the zero-mean 1/f noise signal (pink noise). Frequency range was set to 0.1-10 Hz. Bottom: Representative trace of applied GVS (72-second duration) with an accompanying probability density histogram. Stimulation current magnitude was set to 1 mA and the probability density function follows a Gaussian curve. reported sensing cutaneous stimulation where the electrodes were applied. After determining the stimulus threshold, subjects were seated 80 cm from a screen dis- playing a fixed target, and were instructed to focus their gaze on the target in order to min- imize distractions while the EEG was recorded with and without the delivery of GVS (6 trials). In each trial, EEG was first recorded without the presence of vestibular stimulation for 60 seconds. Noisy stimulation signals were then delivered with a fixed duration of 72 seconds, followed by another 60 seconds of a constant zero current used as a control. Various current intensities were applied in each trial in a pseudo-random manner during the 72-second stimulation period with current levels evenly distributed at fractions of the sensory threshold (10, 26, 42, 58, 74 and 90%). 5.3.4 Data Collection and Preprocessing EEG recordings were conducted using a Neuroscan Synamps2 EEG acquisition system (Neu- roscan Inc., VA, United States). A standard electrode cap (Neuroscan Inc.) was used in con- 76 junction with electrode paste (Electro-Gel; Electrode-Cap International, OH, United States) to decrease overall impedance to less than 10 kΩ in all electrodes. Eighteen EEG channels with one ground and one reference electrode were utilized. A common mastoid reference was used and electrode positions were set according to the standard 10-20 system [128]. Tin cup electrodes (Electrode-Cap International) and surface EMG electrodes were positioned on the ears and eyes, respectively, for artifact removal during analysis. All data were sampled at 1 kHz, with band pass filter settings at 1 to 250 Hz. After artifact removal by ICA using standard methods [124], we examined EEG for both amplitude and connectivity changes. The denoised data were then bandpassed at 1-100 Hz, normalized to unit variance, and de-trended (using the ‘detrend’ function in Matlab). In order to simplify and aid interpretation of the analysis, we averaged the EEG activity in all 20 channels over five ROIs as previously described [109]. The ROIs were grouped as follows: fronto-central (FCentral), left sensorimotor (LSm), right sensorimotor (RSm), central and occipital (see Fig. 3.2 for groupings of electrodes). 5.4 Results 5.4.1 Threshold Test During the threshold test, subjects reported common sensations, such as the feeling of prick- ling or mechanical drumming behind the ear, and did not report hearing any sounds in re- sponse to GVS. Upon becoming aware of the sensory effects at greater current intensities, most individuals began to exhibit slightly greater sensory sensitivity. Therefore, the percep- tual threshold was obtained by decreasing current levels from the high-end threshold value. Eventually, a unique stimulus level was determined at which the subject did not feel the ef- fects of GVS, and this was set generally at a level lower than the subject’s initial high-end perceptual threshold. Subjects reported a sensory threshold at average current amplitudes of 163 ± 110 µA. 77 5.4.2 EEG Analysis Overall, subjects exhibited increased relative amplitudes of EEG activity in the alpha and beta bands in response to rising subthreshold current intensities (Fig. 5.2). The relative power of the beta-band EEG activity increased in response to greater GVS intensities in all ROIs. The correlational trends based on spectral density analysis were significant for relative beta power measured in the LSm, RSm and Occipital regions (p < 0.001, p = 0.004 and p = 0.006, respectively). The relative power of the alpha band modestly increased in response to greater current intensities in the FCentral, LSm, Central and Occipital regions, although the trend was only significant in the FCentral region (p= 0.012). The sparse mAR-PDC based analysis revealed that GVS modulated the strength of direc- tional interaction between pairs of particular ROIs. Connections that were notably enhanced by GVS delivery were: RSm→LSm, LSm→FCentral, RSm→FCentral and Occipital→FCentral (Fig. 5.3). GVS strengthened the RSm→LSm connection with maximal modulation occur- ring right below 20 Hz – that is, in a range comprising alpha and lower beta frequencies. Additionally, connectivity was monotonically augmented in relation to current intensity; yet, this relation did not reach statistical significance. 5.5 Discussion Our results demonstrate that noisy GVS, at an imperceptible level, directly influences brain rhythms in normal, healthy subjects. Furthermore, GVS effects on oscillatory rhythms were frequency-dependent and region-specific. The changes in brain rhythms that we observed are presumably due to direct manipulation of the vestibular system. Galvanic stimulation has been shown to directly alter the firing rates of primary vestibular afferents, thereby bypassing peripheral receptors in the labyrinth [117]. Nevertheless, postural sway at appropriate GVS current intensities still exists in patients who have undergone neurectomy of the vestibular nerve [129], suggesting that GVS may also excite more centrally located vestibular struc- tures. By delivering noisy GVS to healthy subjects, we find an interaction exists between the vestibular system and brain rhythms, which may be mediated by vestibular projections to 78 18 Figure 5.2: Power spectra during application of GVS. Changes in the relative powers of alpha and beta rhythms are plotted as a function of GVS current levels 1-6 (respec- tively, 10, 26, 42, 58, 74 and 90% of sensory threshold) for each ROI (FCentral, LSm, RSm, Central and Occipital). Regression analysis demonstrates that a sig- nificant linear relation exists between the current intensity of GVS and amplitude of beta rhythms in the LSm, RSm and Occipital regions, and the amplitude of alpha rhythms in the FCentral region (p< 0.05). 79 19 Figure 5.3: PDC spectra during application of GVS. Directional connectivity flows from the source of activity represented by each column element towards the sink of activity represented by each row element. Current levels 1 to 6 represent the intensity of galvanic current applied respectively at 10, 26, 42, 58, 74 and 90% of sensory threshold. GVS enhanced the connectivity values of directional flow from LSm, RSm and Occipital regions towards the FCentral region. Connectivity strength for activity from RSm→LSm depended on the frequency of EEG with maximal strength below 20 Hz. 80 multiple thalamic nuclei and/or cortical areas. Numerous animal studies have shown that the vestibular system exerts an influence on subcortical structures, such as the thalamus. Second order afferent neurons not only act as premotor neurons in vestibulo-ocular and vestibulo-spinal reflexes [114], but also project to the cerebellar flocculus and to the thalamocortical pathway for higher-level multisensory cog- nitive processing [117]. The complexity of the vestibular system begins here at the thalamic level: vestibulothalamic projections terminate in multiple thalamic nuclei, where vestibular input is subsequently distributed to widespread cortical areas [116]. For instance, vestibu- lothalamic projections have been shown to be present in the thalamic ventroposterior com- plex, posterior nuclear group (geniculate bodies, pulvinar) and intralaminar thalamic nuclei. Vestibular input from the brainstem nuclei also arrives in the ventroanterior (VA) and ventro- lateral (VL) nuclei, which are connected to the primary motor and premotor cortices [116]. Studies in non-human primates have shown that the VL thalamus receives further input from cerebellar nuclei [130] whereas the VA thalamus receives input from the basal ganglia [131]. The VA-VL complex therefore acts as convergent point where vestibular and motor structures intersect, defining a major vestibulomotor pathway and a conduction course for GVS input. In addition to modulating relayed information, the thalamus functions in basic corticocortical communication; the pulvinar and medial geniculate nuclei contain vestibular-responsive neu- rons, and higher-order relay properties of these thalamic nuclei may allow GVS input to also travel through cortico-thalamo-cortical loops [132]. The numerous projections to thalamic nuclei confer many possible means through which GVS is able to modulate brain rhythms. Our findings of widespread cortical rhythm changes in response to GVS are compatible with previous work in cortical activation. Functional neuroimaging studies in humans have demonstrated that vestibular information is represented in multiple areas: the somatosen- sory cortex, intraparietal sulcus, inferior parietal lobule, posterior insula, temporo-parietal junction, frontal cortex, hippocampus and cingulate cortex [133–136]. While the EEG we used necessarily lacks the spatial resolution of fMRI, we observed similar widespread EEG changes that support the notion that vestibular stimulation affects a distributed network of 81 cortical regions – including those associated with primary sensory, motor and visual function – rather than a restricted region. Our results further demonstrate a pattern of enhanced connectivity directed towards the FCentral region in response to GVS. Kluge et al. [137] showed that epileptic activity in the same area underlies vertigo symptoms, supporting the role of a frontal vestibular processing area. Given that vestibular stimulation activates the premotor and primary motor cortices [134, 135], our results agree with prior studies to demonstrate that vestibular information modulates processing in frontal motor areas. GVS, specifically, has an additional effect to help guide movements [138] in contrast to its effect on stationary postural studies [139]. We note that previous studies have demonstrated that GVS input impacts dynamic motor tasks [140, 141]. Therefore, the beta-band changes in cortical rhythms and connectivity due to GVS deliv- ery may have particular implications for PD. In PD, disruption of abnormal synchronization of the beta-band rhythm in corticostriatal circuitry is associated with therapeutic motor im- provement [142]. For example, high-frequency deep brain stimulation, an advanced surgical treatment for the disease, modulates beta synchrony in the subthalamic nucleus [41]. Tropini et al. [109] found that alterations in beta-band connectivity are a sensitive marker of the dis- ease state of PD patients on and off dopaminergic medication. Our results may provide a pathophysiological mechanism for the previously reported beneficial motor effects of GVS in PD [120, 121]; however, further work will be needed to confirm this. One mechanism through which GVS may affect ongoing cortical rhythms is through “stochastic facilitation”. Stochastic facilitation (previously called “stochastic resonance”) de- scribes the process in which biologically relevant noise may improve information processing in the nervous system [143]. If a non-linear system, such as neuron, is partially depolarized (but not to threshold) by a stimulus, then adding noise to the stimulus may result in random increases and decreases of the stimulus-induced depolarization. Adding noise to the stimu- lus may result in intermittent depolarization as the weak stimulus is now partially detectable [144]. Stochastic facilitation has been suggested as the mechanism through which noisy 82 GVS improves visual memory when non-noisy GVS does not [118]. A similar mechanism may explain our results: noisy stimuli may result in synchronization of neuronal assemblies and phase resetting [144], which may explain the augmentation of (predominantly beta-band) rhythms that we observed. In contrast, modulation of autonomic reflexes has been suggested to be responsible for the improvement in trunk activity and reaction time in a Go/No-Go test in response to GVS [120]. In support of this view, sinusoidal GVS has been shown to influ- ence vasovagal functions, such as blood pressure and heart rate in rats [145]; however, direct autonomic effects of noisy GVS are still unknown. The effects of GVS could potentially be otherwise attributed to activation by top-down mechanisms such as attention and/or general arousal. However, subjects were asked to fo- cus their visual attention on a fixed target in order to minimize distractions. In fact, while attentional mechanisms may modulate vestibular function via corticothalamic feedback, this has only been reported in vestibulo-oculomotor and vestibulo-spinal functions [116]. Addi- tionally, our stimulus was below sensory threshold and resulted in specific effects on brain rhythms; thus, in agreement with previous work [118], we do not believe our observed effects of GVS were a result of either a state of non-specific arousal and/or attentional mechanisms. In summary, we show that imperceptible, noisy GVS specifically affects network oscil- lations, altering both amplitude and directional connectivity in major brain regions in neuro- logically healthy individuals. Moreover, the effects of GVS depend upon the frequency of the oscillations, with notable changes observed in the beta range. Prior studies suggesting a therapeutic effect of GVS on visuo-spatial hemineglect and stroke impairment [146] may be partly mediated by effects on ongoing rhythms. The application of GVS to other neurological conditions known to be affected by abnormal rhythms needs to be more fully explored. 83 Chapter 6 Generalized Multivariate Autoregressive Framework for EEG Source Separation 3 6.1 Introduction Recent advances in neuroimaging techniques including fMRI, MEG and EEG, have made possible the investigation of functional coordination between brain regions during cogni- tive and motor tasks. Brain connectivity may provide a mechanism for functionally related, but spatially distinct neuronal groups to act synergistically. Several mathematical methods have been explored to provide a quantitative measure of brain connectivity using EEG data [9], including coherence, GC [45], PDC [13] and directed transfer function (DTF) [17]. In most studies, these connectivity measures are directly computed from raw scalp EEG signals. However, due to volume conduction, where a deep source of electrical potential can influence many sensors on the scalp surface, the signal measured from an EEG electrode does not ex- clusively represent the activity of one local neural source, but rather the superposition of the activity of several active sources throughout the brain. This can give rise to spurious instan- taneous correlations between scalp EEG signals and potentially lead to misinterpretation of 3A version of this chapter has been published as: J. Chiang, Z. J. Wang, and M. J. McKeown, “A generalized multivariate autoregressive (GmAR)-based approach for EEG source connectivity analysis,” IEEE Transactions on Signal Processing, vol. 60, no. 1, pp. 453–465, Jan. 2012. c© 2012 IEEE. http://dx.doi.org/10.1109/TSP.2011.2166392 84 the connectivity results [75]. To address the problem of volume conduction in EEG, several new connectivity analysis techniques have been proposed. These include computing the surface Laplacian as a spatial “whitening” filter to decorrelate activity in neighboring channels [147], and, more recently, studying the brain connectivity in the “source space” [7, 148, 149]. The source space meth- ods attempt to estimate the underlying neural sources from scalp EEG signals through the inversion of the mixing process of volume conduction and subsequently apply the standard connectivity measures to the source estimates. In general, determining neural sources from EEG signals is a mathematically ill-conditioned inverse problem, and it has no unique solution without prior knowledge or statistical assump- tions [150]. One common approach to solving the source identification problem is the dipole modeling technique, which aims to estimate the electrical dipoles that best describe the ob- served EEG signals based on the assumed forward head model. The successful estimation of the dipole sources requires proper specification of the head model including the number of ac- tive dipoles, models of the generators of the electrical activity, and the shape and conductivity of the head model [22]. Another possible approach to finding sources in the EEG is to use blind source separation (BSS) techniques [151–153]. BSS aims to recover the source signals, as well as the unmixing transformation, from the observed mixed signals by exploiting a priori knowledge about the statistical properties or structures of the sources. BSS techniques can be appealing as it separates the problem of EEG source identification from the problem of source localization which requires an explicit specification of forward head model. One of the most widely used BSS algorithms is ICA, which assumes the mutual independence of underlying sources [23]. It has been commonly used to isolate eye movement and muscle noise artifacts in EEG [124], as these artifacts are typically independent from ongoing EEG brain rhythms. However, the underlying assumption of statistical independence between the activations of different neural assemblies still remains to be validated experimentally [150]. Furthermore, as we try to establish the causal relationships and connectivity between EEG sources, it is 85 implicitly assumed that the sources may be influenced by each other. This contradicts the fundamental assumption of mutual independence between sources in ICA. Recently, Gómez-Herrero et al. propose an integrated approach for estimating the source connectivity from scalp EEG measurement by formulating the volume conduction effects and the connectivity between underlying brain sources as a state-space mAR framework [148]. In this framework, the sensor measurement is assumed to be generated from an instantaneous linear mixing of brain source activity whose spatial and temporal structure is modeled by an mAR process. With the assumption of non-Gaussian innovation for the mAR process at the source level, a two-step estimation procedure, named MVAR-ICA, is developed which determines the unknown mixing matrix and mAR parameters via the sequential applications of mAR fitting to EEG signals and ICA demixing to the residuals of the fitted mAR. However, the performance of this two-step approach is highly dependent on the proper estimation of the mAR model as any error incurred in this step gets carried over to the later steps and affects the final estimation of source connectivity. Other related works extending the general state-space mAR framework with additional model assumptions or parameter constraints are developed. In [149], the authors incorporate sparsity constraints on mAR coefficients and propose a convolutive ICA-based estimation technique. In [154], a cortical patch basis representation is used to describe the source activity distribution within each predefined ROI at the sensor level and the basis coefficients for each ROI are estimated using the expectation maximization algorithm. In this study, we extend the state-space model proposed in [148] and formulate the causal relationships between underlying brain source as a generalized mAR process with unknown AR coefficients and noise parameters. Here, the term “generalized mAR” means that the innovation process of the mAR model follows a generalized Gaussian distribution, which al- lows more flexible modeling of sub- and super-Gaussian signals, as those commonly observed in neurophysiological signals [155]. In contrast to the conventional two-stage approaches, we propose a ML procedure to simultaneously estimate both the mixing matrix of volume con- duction and mAR model parameters directly from EEG signals. We refer to this estimation 86 technique as “Generalized mAR Source Separation” or GmAR-SS. In line with the aforemen- tioned works [148, 149, 154], we assume that all information flows between brain sources are time-lagged (at plausible sampling rates) due to propagation delays and there is no simulta- neous source co-activation. We note that with these assumptions it is possible to distinguish between source interactions of interest from the potential serious confound of volume con- duction. While these assumptions will need ultimate verification, they are, however, no more stringent or unrealistic than those of, e.g., ICA. The connectivity estimation performance of the proposed GmAR-SS technique is verified with simulated EEG data and compared with existing two-stage approaches. The proposed approach is also applied to real EEG data collected from 8 healthy subjects and 9 PD patients during a right-handed bulb-squeezing task with differing amounts of visual feedback. We demonstrate that the proposed approach is able to extract robust sources whose connectivity differs across groups and task conditions. 6.2 Methods In this section, we first describe the state-space generalized mAR model and the ML esti- mator. Next, we present the group analysis technique for real EEG data, which is adapted from the group analysis procedure proposed in [148]. Finally, we describe the simulation framework for comparing performance of various connectivity estimation techniques. 6.2.1 EEG Signal Model First, we start with the state-space mAR framework proposed in [148]. Let x(t) = [ x1(t) x2(t) . . . xJ(t) ]T denote an J-dimensional EEG data vector at time t, where the superscript T denotes matrix transpose. We assume that x(t) is generated from an instantaneous linear mixing of the sources s(t), which can be written as x(t) =Cs(t), 1≤ t ≤ N (6.1) 87 where s(t) = [ s1(t) s2(t) . . . sM(t) ]T denotes the underlying EEG sources, and C is an J×M mixing matrix whose i-th column represents the projection from the i-th source to the scalp electrodes, which is also commonly referred to as “scalp projection map” in the ICA literature. Note that, here, we focus on the case where the number of sources is equal to the number of sensors, i.e., J =M. At the source level, we model s(t) by a generalized mAR process of order P, which is written as s(t) = P ∑ p=1 Aps(t− p)+v(t), (6.2) where Ap ∈ RM×M, p = 1, . . . ,P are the mAR coefficient matrices capturing the causal rela- tionships between brain sources, and v(t) = [ v1(t) v2(t) . . . vM(t) ]T denotes the inno- vation process whose elements are assumed to be mutually independent. Different from the non-Gaussian assumption made in [148, 149], we assume that each element of the innovation process follows a zero-mean, univariate generalized Gaussian distribution with the density p(vi(t)|Ri,wi) = Ri2wiΓ(1/Ri)e −| vi(t)wi | Ri , (6.3) where Γ(·) is the gamma function and the noise parameters Ri and wi determine the shape (peakiness) and the width of the distribution, respectively. The generalized Gaussian distri- bution is a more general parametric form of a Gaussian which can model a super-Gaussian distribution with Ri< 2 (e.g., Laplacian), a Gaussian with Ri= 2, and a sub-Gaussian distribu- tion with Ri> 2 (e.g., uniform). Because of it flexibility, the generalized Gaussian distribution has been used in several different ICA models for modeling EEG signals [155, 156]. In summary, the proposed state-space generalized mAR model can be specified by the following parameter set θ = {C,Ap,Ri,wi}, 1≤ p≤ P, 1≤ i≤M. (6.4) 88 Our objective is to, given the observed EEG signals x(t), estimate the model parameter set θ , together with the time courses of s(t). 6.2.2 Maximum Likelihood Estimation The model parameter set θ is estimated using an iterative ML approach. The basic idea is that, at each iteration, the unobserved sources s(t) are first computed using the unmixing matrix W (defined below) estimated from the previous iteration, and the sources are then used, as if they were observed, to update the log-likelihood function. The model parameters are re-estimated based on the updated log-likelihood function. The iteration repeats until convergence is reached. Let us first define the data likelihood (likelihood of observations) with respect to θ as L, px(x(1),x(2), . . . ,x(N)|θ ). (6.5) LetW denote the unmixing matrix which is defined asW =C−1 such that s(t)=Wx(t). Using the multivariate transformation [157], the likelihood of observations defined in Eqn. 6.5 can be calculated from the likelihood of the source densities: L = ( 1 |det(W−1)|) N ps(Wx(1),Wx(2), . . . ,Wx(N)) = (|det(W )|)N ps(s(1),s(2), . . . ,s(N)), (6.6) where | · | denotes the absolute value. Since the sources are assumed to follow a P-th order mAR process, the conditional prob- ability of each time point depends only upon its previous P time points, i.e., ps(s(t)|s(t−1), . . . ,s(1)) = ps(s(t)|s(t−1), . . . ,s(t−P)) = pv(s(t)− P ∑ p=1 Aps(t− p)). (6.7) 89 Under the mutual independence assumption of the residuals v(t), Eqn. 6.7 can be simplified as ps(s(t)|s(t−1), . . . ,s(1)) = M ∏ i=1 pvi(e T i (s(t)− P ∑ p=1 Aps(t− p))|wi,Ri), (6.8) where ei is the i-th column of an M×M identity matrix. Combining Eqn. 6.6 and Eqn. 6.8, we can rewrite the data likelihood as follows L= (|det(W )|)N ps(s(1), . . . ,s(P)) (6.9) · N ∏ t=P+1 M ∏ i=1 pvi(e T i (s(t)− P ∑ p=1 Aps(t− p))|wi,Ri). By taking the logarithm of Eqn. 6.9, we obtain the log likelihood function of the observations LL , log px(x(1),x(1), . . . ,x(N)) = N log(|det(W )|)+ log ps(s(1), . . . ,s(P)) + N ∑ t=P+1 M ∑ i=1 [ log Ri 2wiΓ(1/Ri) − ∣∣∣∣zi(t)wi ∣∣∣∣Ri ] , (6.10) where zi(t) = eTi (s(t)− P ∑ p=1 Aps(t− p)). (6.11) The ML estimates of model parameters are found by setting the respective partial deriva- 90 tive of the log-likelihood (i.e., ∂LL∂W , ∂LL ∂Ap , ∂LL ∂Ri and ∂LL∂wi ) to zero. This gives wi = ( Ri N N ∑ t=P+1 |zi(t)|Ri ) 1 Ri , (6.12) ∂LL ∂Ri = N Ri + N R2i Γ′( 1Ri ) Γ( 1Ri ) − N ∑ t=P+1 [( |zi(t)| wi )Ri log |zi(t)| wi ] = 0, (6.13) ∂LL ∂Ap(i, j) = Ri wRii N ∑ t=P+1 |zi(t)|Ri−1sign(zi(t))s j(t− p) = 0, (6.14) ∂LL ∂Wi, j = NWi, j− Ri wRii |zi(t)|Ri−1sign(zi(t))eTjW Ts(t) = 0. (6.15) Since Eqns. 6.13, 6.14 and 6.15 do not have closed-form solutions, we solve them numer- ically using a gradient ascent approach such as the Newton-Raphson method. It is important to note that as the learning is based on the ML framework, the algorithm may only converge to a local maximum rather than a global maximum. To address this issue, we utilize a multi- start strategy where multiple random initial conditions are used to generate multiple solutions and the one with highest likelihood is chosen as the final estimate. 6.2.3 Group Analysis An important concern in biomedical research is group analysis, since it is desirable to extract sources from subjects that are representative of the group as a whole. However, when the source separation is applied to each subject separately, the ambiguity in the ordering of the extracted sources makes the distillation of results across subjects difficult. To address this issue, we adopted the group analysis approach presented in [148, 158] by performing source separation on each subject separately and using the clustering technique from Icasso [159] to find source components that have high repeatability across subjects. An illustration of the group analysis for finding group-level connectivity pattern is shown in Fig. 6.1. The group- level estimation works by clustering all subjects’ estimated sources (including both normal and PD subjects) given by GmAR-SS according to the correlation coefficients between their scalp projection patterns (i.e., columns ofCelec) using agglomerative clustering. The number 91 Subject-level Estimation x(t)Subj 1 x(t)Subj K . . . . GmAR-SS and and AP A1 AP A1 … Group-level Estimation Similarity Matrix Clustering Φ … … … Connectivity Estimation Pseudo- inverse Ψcommon s*(t)Subj 1 s*(t)Subj K . . . . Multiply x(t)Subj 1 PDC GmAR-SS … Multiply x(t)Subj K 10 20 30 40 50 60 70 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Frequency (Hz) P D C Frequency EEG signals of subjects AR coefficient matrices Scalp projection maps Similarity between scalp projection maps Centrotypes of all clusters Group spatial filter Source signals of subjects Connectivity estimates Select rows corr. to significant clusters Figure 6.1: Estimation procedure for group connectivity analysis. This figure is prepared based on the group analysis framework originally proposed in [148]. 92 of clusters in the data is determined visually based on the size and tightness of clusters. As some clusters may be associated with spurious activity that only appears in a single or small number of subjects, for subsequent connectivity analysis, only significant clusters that have high inter-subject repeatability are considered. Specifically, we consider a cluster as significant if it contains sources from at least 50% of the subjects. Each of the identified cluster is represented by a “centrotype”, which is the source estimate in the cluster that has the maximum sum of correlation coefficients with the rest of the cluster members. Assuming that the brain sources activated by the performed task and their scalp projection patterns do not vary considerably across subjects, the source ambiguity issue inherent in source separation techniques can be solved by deriving a common, global spatial filter [148], which allows to extract source activities that are comparable across subjects. Let Φ denote the group mixing matrix formed by concatenating the centrotypes of all clusters: Φ = [φ 1, . . . ,φQ], (6.16) where Q is the total number of clusters and φ q ∈R19×1 denotes the scalp projection pattern of the q-th centrotype. The group spatial filter, denoted byΨcommon, is constructed by selecting rows of the pseudo-inverse ofΦ that correspond to the significant clusters. Subsequently, for any arbitrary subject, the temporal dynamics of the brain sources can then be recovered from the EEG recordings by s∗(t) =Ψcommonx(t), (6.17) where x(t) denotes the EEG measurements and s∗(t) denotes the sources recovered by spatial filtering. 6.2.4 Simulations The estimation performance of the proposed GmAR-SS technique is evaluated using sim- ulated EEG data generated based on a single-shell spherical head model in [160]. Four (M = 4) source dipoles are simulated whose activity, s(t), follows a 5th-order (P = 5) gen- 93 drm θm z Right Left x y βm em Dm Figure 6.2: Spatial characteristics of a source dipole. The dipole is specified by three parameters: (1) radius rm, which is the distance between the dipole and the center of the head, (2) spatial range parameter, βm, which is the angle between the radial component of the dipole and the z-axis, and (3) orientation angle, θm, which is the angle between the dipole and its radial component. eralized mAR process. The mAR coefficient matrix Ap are randomly generated in such a way that the resulting mAR process is stable and the ratio of non-zero off-diagonal entries in Ap,∀p = 1, . . . ,P, to the total number of off-diagonal entries is specified by c. The inno- vations vm(t), m = 1, . . . ,M, driving the generalized mAR process are mutually independent and each of them is drawn from a generalized Gaussian distribution whose probability den- sity function is specified by the shape parameter Rm and scale parameter wm. The position and orientation of each source dipole are specified by three parameters (see Fig. 6.2): (1) the radius, rm ∈ [0,1], which is the distance between the dipole and the center of the head (normalized to head radius), (2) the spatial range parameter, βm ∈ [−45◦,45◦], which is the angle between the radial component (−→e m in Fig. 6.2) of the dipole and the z-axis, and (3) the orientation angle, θm ∈ [0◦,90◦], which is the angle the dipole −→D m makes with its radial component −→e m (e.g., θm = 0◦ if −→D m is perpendicular to the scalp and θm = 90◦ if tangential to the scalp). Sixteen EEG electrodes (J = 16) are simulated which are evenly spaced on the surface of the simulated head with normalized radius of 1 (see the scalp projection maps in Fig. 6.3). Voltage gains (i.e., scalp projections) of the source dipoles at the simulated elec- trodes are determined by Brody’s formula [161] and they form columns of the mixing matrix 94 C. The simulated EEG signals, x(t), are generated by multiplying C by the source dipole activity. To evaluate the robustness of the proposed technique to the presence of unmodeled noise, in all of our simulations, we also add 4 independent noise dipoles whose individual activity, un(t),n = 1, . . . ,4, follows a 3rd-order univariate AR model driven by white Gaussian noise. These noise dipoles are added to represent the biological noise (possibly coming from the brain or elsewhere, e.g., heart) which are unconnected to the interacting brain sources of interest. The positions and orientations of the noise dipoles, specified by the parameters rnoisen , β noisen and θ noisen , are randomly varied in each simulation with values drawn uniformly from the following ranges: 0.3 ≤ rnoisen ≤ 0.9, −45◦ ≤ β noisen ≤ 45◦ and 0◦ ≤ θ noisen ≤ 90◦. The voltage gains of the noise dipoles at simulated electrodes are represented by the noise mixing matrix, Cnoise. We also simulate sensor noise, n(t), by adding white Gaussian noise to the simulated EEG signals. In summary, the signals measured at the electrodes can be expressed by x̃(t) =Cs(t)+Cnoiseu(t)+n(t), (6.18) where u(t) = [u1(t), . . . ,u4(t)]T is the noise dipole activity. The variance of the biological noise is determined based on the desired signal-to-biological noise ratio (SBNR), which is defined as the average standard deviation of clean EEG signalsCs(t) divided by the standard deviation of the biological noise Cnoiseu(t). Similarly, the variance of sensor noise is deter- mined by the signal-to-sensor noise ratio (SNR), which is defined as the average standard deviation ofCs(t) divided by the standard deviation of n(t). To test the estimation performance and reliability of the proposed method in identify- ing source connectivity, we conduct 4 simulation experiments to investigate the effects of the following factors: 1) SNR, 2) SBNR, 3) sample length (N) and 4) R value (Gaussian- ity) of the source innovations. In each simulation experiment, 50 sets of EEG signals, each with a different realization of AR coefficients, dipole locations and orientations, source innovations and noise, are generated for every input value of the tested parameter. The 95 parameters that are not being tested are set either to some default values (SNR = 10dB, SBNR = 10dB, N = 10,000, c = 10/12) or to values drawn randomly from some plausi- ble ranges (1≤ Rm ≤ 3, 0.1≤ wm ≤ 1, 0.3≤ rm ≤ 0.9, −45◦ ≤ βm ≤ 45◦ and 0◦ ≤ θm ≤ 90◦ for all m = 1, . . . ,M) unless specified otherwise. Such arbitrary combinations of parameter values enable us to obtain a statistically meaningful evaluation of the estimation performance of the proposed technique under different scenarios. During the estimation procedure, the mAR model order is assumed unknown and the automatic order selection technique, Akaike information criterion (AIC), is used to determine the optimal mAR order from the range of 1 to 20 in all simulations. The resulting simulated EEG data x̃ is of dimension J×N. Prior to applying the GmAR- SS algorithm, we reduce the dimension of x̃ to M×N using a combination of principal compo- nent analysis (PCA) and partial least square (PLS) [162]. PCA is one of the most commonly used dimension reduction techniques which projects high-dimensional data onto a subset of principal component directions that capture the most variance in the data. However, our pre- liminary study showed that using PCA alone may potentially remove useful information in the data and leads to degradation in estimation performance. Thus, we propose combining PCA with PLS in which we model the simulated EEG data as predictors and the entire set of the principal components as responses. The goal of PLS is to iteratively extract orthog- onal linear combinations of predictors that best explain variance in the predictor data and, at the same time, have maximal correlations with the response data [162]. The first M PLS factors are then analyzed by the proposed GmAR-SS algorithm. Finally, the original mixing matrix is recovered by multiplying the inverse of the PLS loading matrix with the estimated unmixing matrix. Performance measures To measure the performance of unmixing (i.e., the accuracy of C estimates), the following metric [155] is used Eunmixing = ||ŴCP−diag(ŴCP)|| ||ŴC|| , (6.19) 96 whereC denotes the true mixing matrix and Ŵ is the estimated unmixing matrix. The matrix P is a permutation matrix which is included to account for the ambiguity in the ordering of the estimated sources and it is determined using the pairing technique proposed in [163]. To evaluate the estimation accuracy of the mAR coefficient matrix Ap, we compare the PDC obtained from the true Ap with those given by the tested method. The PDC technique is chosen because it provides normalized values bounded between 0 and 1, which simplifies the evaluation of Âp. The estimation errors of the source PDC are measured by the following metric [148]: EPDC = 1 M2F ∑i, j∈{1,...,M} f=0:0.1:0.5 √ (pi j←i( f )− p̂i j←i( f ))2, (6.20) where pi j←i( f ) denotes the PDC value given by the true Ap, and p̂i j←i( f ) is the PDC given by the mAR coefficient estimate Âp. F is the total number of frequency bins used in the PDC calculations (in this study, F = 51). Comparative schemes The performance of the proposed GmAR source separation technique is compared with three other techniques: the traditional mAR approach, ICA and MVAR-ICA. The traditional mAR approach fits an mAR model directly to the observed signals x(t), assuming that they are identical to the underlying sources. This essentially disregards the mixing effects caused by volume conduction. To estimate the mAR parameters, the ARfit algorithm developed in [164] is used. The second tested technique is the ICA approach which involves first applying ICA to EEG data to remove the volume conduction effects and subsequently fitting an mAR model to the extracted independent components to infer source connectivity. The third tested technique, MVAR-ICA method [148], models measured EEG signals as an instantaneous linear mixture of unobserved source signals. The sources are assumed to follow an mAR process whose innovations are modeled by mutually independent non- Gaussian random variables. Assuming that there is no noise at the sensor level, the observed 97 EEG x(t) are rewritten in the form of an mAR process: x(t) = P ∑ p=1 CApC−1x(t− p)+Cv(t) (6.21) = P ∑ p=1 Bpx(t− p)+Cv(t), (6.22) where C denotes an unknown mixing matrix modeling volume conduction from the sources to the scalp electrodes and v(t) denotes the innovation process of the source-level mAR. The matrices Ap measure the causal relationships between sources whereas B̂p ,CApC−1 measure the relationships between EEG signals. MVAR-ICA solves the identification problem by first fitting an mAR model to the EEG signals using ARfit [164]. The estimated AR matrices B̂p capture the time-lagged flows of information that is produced by propagation of neural activity. The remaining instantaneous cross-dependencies are assumed to be exclusively caused by volume conduction and are cap- tured in the residuals of the fitted mAR model, Cv(t). Given that the elements of v(t) are mutually independent, the mixing matrixC can be recovered by applying ICA toCv(t). The estimated ICA mixing matrix Ĉ is used to estimate the mAR coefficient matrices as: Âp = Ĉ −1 B̂pĈ. (6.23) Note that prior to application of the aforementioned comparative schemes, dimension reduc- tion is first performed using PCA [148] and only the first M principal components are kept for subsequent source separation analyses. 98 0 0.2 0.4 0 0.5 1 S2 −> S1 0 0.2 0.4 0 0.5 1 S3 −> S1 0 0.2 0.4 0 0.5 1 S4 −> S1 0 0.2 0.4 0 0.5 1 S1 −> S2 0 0.2 0.4 0 0.5 1 S3 −> S2 0 0.2 0.4 0 0.5 1 S4 −> S2 0 0.2 0.4 0 0.5 1 S1 −> S3 0 0.2 0.4 0 0.5 1 S2 −> S3 0 0.2 0.4 0 0.5 1 S4 −> S3 0 0.2 0.4 0 0.5 1 S1 −> S4 0 0.2 0.4 0 0.5 1 S2 −> S4 0 0.2 0.4 0 0.5 1 S3 −> S4 Freq P D C True PDC ICA MVAR−ICA GmAR−SS Figure 6.3: PDC estimation results from an example simulation run. Diagonal subfigures: Scalp projection maps of simulated sources. Off-diagonal subfigures: PDC spectra estimated using different source separation techniques. True PDC spectra (with crosses) are also shown. 99 (a) S1 (b) S2 (c) S3 (d) S4 Figure 6.4: Scalp projection maps of the estimated sources given by GmAR-SS from the same simulation run as Fig. 6.3. The sign and order of the estimated sources are permuted to match maximally with the projection maps of simulated sources. 6.3 Simulation Results For all of the simulations results presented in this section, default parameter values/ranges described in Section 6.2.4 are used unless specified otherwise. In addition, the simulated EEG signals are reduced to exactly as many dimensions as the simulated sources (i.e., J =M = 4) before the connectivity analysis techniques are applied. The results of an example simulation run are shown in Fig. 6.3 and Fig. 6.4. In this example, the level of connectedness of the source mAR coefficient matrix, c, is set to 33% (i.e., 4/12 of all possible interactions presents). The scalp projection patterns describing the influence of simulated source dipoles on the EEG (i.e., columns of C) are shown in the diagonal subfigures of Fig. 6.3. The scalp projections of the reconstructed sources given by the proposed GmAR-SS technique (after the sign and order are permuted to match the simulated sources) are shown in Fig. 6.4. It can be seen that GmAR-SS achieves great unmixing matrix estimation performance, yielding an unmixing error at 18.8%. In regards to the PDC estimation (see Fig. 6.3), GmAR-SS achieves the lowest PDC error at 10.5%, compared to the traditional mAR at 12.1% (not shown), ICA at 14.1% and MVAR-ICA at 16.9%. The performance improvement of GmAR-SS is most noticeable in connections S1→ S2 and S1→ S4, for which among all the tested methods, only GmAR-SS is able to correctly identify both the location and strength of the spectral peaks in the PDC spectra. Fig. 6.5 shows the robustness of the tested methods to additive white Gaussian noise 100 at the electrodes. The SNR is varied from 0 dB to 20 dB and the results are obtained by averaging over 50 different realizations. We note that the traditional mAR approach generally yields the highest PDC errors among the tested methods as it ignores the spurious correlations between electrodes caused by volume conduction. For the ICA approach, due to the presence of temporal correlations in the underlying sources, the independence assumption of ICA is violated, leading to poor estimations of the unmixing matrix and mAR matrices. On the other hand, by taking into account the temporal structure of the sources, both MVAR-ICA and GmAR-SS show significant improvements over the ICA approach. In particular, the proposed GmAR-SS generally achieves the lowest PDC errors despite the presence of unmodeled noise at the sensor level. Fig. 6.6 shows the effect of biological noise on mixing matrix and PDC estimation. The relative degradations in the estimation performance of the tested methods with decreasing SBNR show similar trends as in the varying SNR. One noticeable difference is the PDC estimation performance of the GmAR-SS method at low SBNR. When the SBNR is dropped below 2 dB, the PDC error of GmAR-SS increases rapidly compared to MVAR-ICA which is comparably less sensitive to varying SBNR. It is interesting to note that the PDC errors of tradition mAR, ICA and MVAR-ICA ap- proaches are relatively insensitive to the varying noise level (see Figs. 6.5b and 6.6b). A possible explanation for this may be the use of PCA for dimension reduction as the first step of the connectivity estimation. PCA works most effectively when the data are Gaus- sian distributed. However, given that the underlying sources are simulated using generalized Gaussian distribution, this mismatch between the PCA assumption and the simulated source model can potentially cause useful information being partially eliminated during the dimen- sion reduction step. The effect of the sample length (N) on estimation performance of the tested methods is illustrated in Fig. 6.7. In general, the estimation performances of all methods improve as the number of samples increases. In particular, the proposed GmAR-SS achieves the best PDC estimation performance for all considered sample length, even with only 1000 data points. 101 We also note that the performance of MVAR-ICA degrades considerably when the sample length drops below 5000 data points. To investigate the effect of Gaussianity of the mAR innovations, we simulated sources with Rm varying in the range from 1 to 3. Results in Fig. 6.8 show that when the source inno- vations approach a Gaussian distribution (i.e., Rm = 2), the estimation performance of ICA, MVAR-ICA and GmAR-SS deteriorates drastically due to the violation of the assumption of non-Gaussian innovations. As Rm gets further away from 2, the performance of MVAR- ICA and GmAR-SS becomes considerably better. ICA, however, suffers greatly when sub- Gaussian innovations (Rm > 2) are used. The traditional mAR estimator is comparably less sensitive to the distributions of source innovations. In most simulation scenarios, the proposed GmAR-SS and MVAR-ICA achieved com- parable performance in estimating the mixing matrices. However, in estimating Ap, MVAR- ICA yielded a noticeably larger error in PDC relative to GmAR-SS. This intuitively can be explained as follows: MVAR-ICA estimates the source-level mAR coefficient matricesAp us- ing a multi-step approach, which involves sequential estimations of sensor-level mAR model, mixing matrix and then finally source-level mAR model (see Section 6.2.4). Any errors in- curred in the early steps can be propagated into the final estimation of Ap, and thus results in larger PDC error than GmAR-SS. In the simulations, data reduction was done prior to fitting the model and the dimensional- ity of the underlying sources were assumed known a priori. However, in real applications, the source dimensionality is usually unknown. Thus, we next explore the effect of a dimension- ality mismatch between the simulated source space and the dimension-reduced sensor space. To show the effect of a dimensionality mismatch in the dimension reduction procedure, we have taken the same simulation run as the one used to generate Fig. 6.3 and performed PCA to reduce the 16-channel simulated EEG data to 6 dimensions (while the true number should be 4). Our proposed GmAR-SS algorithm was applied to the dimension reduced data and the results are shown in Fig. 6.9 and Fig. 6.10. We can see from Fig. 6.9 that among the 6 estimated sources, S1 and S6 have very simi- 102 0 5 10 15 20 0 10 20 30 40 50 60 70 80 90 100 SNR (dB) Un m ixi ng E rro r ( %) Unmixing Error ICA MVAR−ICA GmAR−SS (a) Unmixing Error 0 5 10 15 20 0 2 4 6 8 10 12 14 16 18 SNR (dB) PD C Er ro r ( %) PDC Error Traditional mAR ICA MVAR−ICA GmAR−SS (b) PDC Error Figure 6.5: Effects of sensor noise, n(t), on mixing matrix and PDC estimation. The sensor noise is simulated by adding white Gaussian noise to each electrode. Bio- logical noise is also present with a fixed SBNR of 20 dB. (a) Average unmixing error Eunmixing. (b) Average PDC error EPDC. 103 0 5 10 15 20 0 10 20 30 40 50 60 70 80 90 100 SBNR (dB) Un m ixi ng E rro r ( %) Unmixing Error ICA MVAR−ICA GmAR−SS (a) Unmixing Error 0 5 10 15 20 0 2 4 6 8 10 12 14 16 18 SBNR (dB) PD C Er ro r ( %) PDC Error Traditional mAR ICA MVAR−ICA GmAR−SS (b) PDC Error Figure 6.6: Effects of biological noise, Cnoiseu(t), on mixing matrix and PDC estima- tion. The biological noise is simulated using 4 randomly placed noise dipoles whose individual activity is generated by a 3rd-order univariate AR process. Sen- sor noise is also present with a fixed SNR of 10 dB. (a) Average unmixing error Eunmixing. (b) Average PDC error EPDC. 104 1 2 3 4 5 6 x 104 0 10 20 30 40 50 60 70 80 90 100 Sample Length Un m ixi ng E rro r ( %) Unmixing Error ICA MVAR−ICA GmAR−SS (a) Unmixing Error 1 2 3 4 5 6 x 104 0 2 4 6 8 10 12 14 16 18 Sample Length PD C Er ro r ( %) PDC Error Traditional mAR ICA MVAR−ICA GmAR−SS (b) PDC Error Figure 6.7: Effects of sample length, N, on mixing matrix and PDC estimation. (a) Average unmixing error Eunmixing. (b) Average PDC error EPDC. 105 1 1.5 2 2.5 3 0 10 20 30 40 50 60 70 80 90 100 R (Gaussianity of Source Innovations) Un m ixi ng E rro r ( %) Unmixing Error ICA MVAR−ICA GmAR−SS (a) Unmixing Error 1 1.5 2 2.5 3 0 2 4 6 8 10 12 14 16 18 R (Gaussianity of Source Innovations) PD C Er ro r ( %) PDC Error Traditional mAR ICA MVAR−ICA GmAR−SS (b) PDC Error Figure 6.8: Effects of Gaussianity of source innovations, Rm,∀i = 1, . . . ,4, on mixing matrix and PDC estimation. (a) Average unmixing error Eunmixing. (b) Average PDC error EPDC. Note the poor performances of ICA, MVAR-ICA and GmAR- SS when the source innovations approach to a Gaussian distribution (Rm = 2). 106 (a) S1 (b) S2 (c) S3 (d) S4 (e) S5 (f) S6 Figure 6.9: Scalp projection maps of the 6 source dipoles recovered by GmAR-SS while the signals were simulated from 4 underlying sources. lar scalp projection maps and S5 has negligible contributions to the activity measured at the electrodes. Thus, upon visual inspection, we would expect that the simulated EEG data were actually generated from 4 source dipoles. The PDC spectra corresponding to the estimated source dipoles 4, 1, 3 and 2 (note that the order of estimated sources is permuted to match that of true sources) are shown in Figure 3 against the true PDC spectra. In this example, GmAR-SS provides a close approximation of the true PDC spectra despite the dimensional- ity mismatch during the dimension reduction step. In addition, to determine the connectivity patterns between the estimated sources, we have constructed 1% significance threshold using surrogate data. The surrogate data were generated by independently and randomly shuffling the temporal order of the estimated source signals. An mAR model was fitted to the surrogate data and the corresponding PDC was calculated. The PDC values from 100 sets of surrogate data were used to construct an empirical distribution from which the 1% significance level is derived (see Fig. 6.11 below). The connectivity patterns between estimated sources and the true sources are shown in Fig. 6.12. We can see that the estimated source S5 is isolated indicating that it is a noisy source. Also, given that the estimated sources S1 and S6 are 107 probably identical sources, they share same connectivity structure with the rest of the esti- mated sources. If we map the estimated sources S1, S2, S3 and S4 to the simulated sources S2, S4, S3 and S1 respectively, we obtain almost the same connectivity pattern except for the connection from estimated source S1 to S4 which is incorrectly identified. Nonetheless, if we compare the PDC patterns in Fig. 6.10 (based on 6 estimated sources) with Fig. 6.3 (based on 4 estimated sources), we can still notice a slight increase in PDC error caused by dimensionality mismatch. 6.4 Application to Real EEG Data - A Case Study 6.4.1 Experimental Procedures and Data Acquisition The proposed method was applied to real EEG data recorded from normal subjects and sub- jects with Parkinson’s disease performing motor tasks. The study was approved by the Uni- versity of British Columbia Ethics Board and all subjects gave informed written consent prior to participating. Nine PD patients (mean age: 66 yrs) displaying mild to moderate level of PD (stage 1-2 according to Hoehn and Yahr assessment) were recruited from Pacific Parkin- son’s Research Centre (Vancouver, Canada). Their motor symptoms were assessed using the UPDRS, with a mean score of 23. Eight age-matched healthy subjects were recruited as controls. During the experiment, subjects were seated 2 m away from a large computer screen. Vi- sual target was displayed on the screen as a vertical yellow bar oscillating in height at one of the following three frequencies: 0.2, 0.4 or 0.6 Hz. Subjects were asked to squeeze a pressure- responsive bulb with their right hand such that the force produced by the subject followed the height of target bar as close as possible. Visual feedback representing the force output of the subject was displayed as a vertical green bar superimposed on the target bar (see Fig. 6.13). The squeezing task was performed under two different conditions: an external guidance con- dition where both the visual target and visual feedback of the force generated were shown (follow); and a self-guided condition where visual feedback was eliminated but the visual tar- 108 0 0.2 0.4 0 0.5 1 S2 −> S1 0 0.2 0.4 0 0.5 1 S3 −> S1 0 0.2 0.4 0 0.5 1 S4 −> S1 0 0.2 0.4 0 0.5 1 S1 −> S2 0 0.2 0.4 0 0.5 1 S3 −> S2 0 0.2 0.4 0 0.5 1 S4 −> S2 0 0.2 0.4 0 0.5 1 S1 −> S3 0 0.2 0.4 0 0.5 1 S2 −> S3 0 0.2 0.4 0 0.5 1 S4 −> S3 0 0.2 0.4 0 0.5 1 S1 −> S4 0 0.2 0.4 0 0.5 1 S2 −> S4 0 0.2 0.4 0 0.5 1 S3 −> S4 Freq P D C True PDC GmAR−SS Figure 6.10: PDC of the selected estimated sources (blue line with circle) and simulated sources (green line with cross). Note that the sign and order of the estimated sources are permuted to match maximally with the projection maps of simulated sources. 109 0 0.5 0 0.5 1 S2 −> S1 0 0.5 0 0.5 1 S3 −> S1 0 0.5 0 0.5 1 S4 −> S1 0 0.5 0 0.5 1 S5 −> S1 0 0.5 0 0.5 1 S6 −> S1 0 0.5 0 0.5 1 S1 −> S2 0 0.5 0 0.5 1 S3 −> S2 0 0.5 0 0.5 1 S4 −> S2 0 0.5 0 0.5 1 S5 −> S2 0 0.5 0 0.5 1 S6 −> S2 0 0.5 0 0.5 1 S1 −> S3 0 0.5 0 0.5 1 S2 −> S3 0 0.5 0 0.5 1 S4 −> S3 0 0.5 0 0.5 1 S5 −> S3 0 0.5 0 0.5 1 S6 −> S3 0 0.5 0 0.5 1 S1 −> S4 0 0.5 0 0.5 1 S2 −> S4 0 0.5 0 0.5 1 S3 −> S4 0 0.5 0 0.5 1 S5 −> S4 0 0.5 0 0.5 1 S6 −> S4 0 0.5 0 0.5 1 S1 −> S5 0 0.5 0 0.5 1 S2 −> S5 0 0.5 0 0.5 1 S3 −> S5 0 0.5 0 0.5 1 S4 −> S5 0 0.5 0 0.5 1 S6 −> S5 0 0.5 0 0.5 1 S1 −> S6 0 0.5 0 0.5 1 S2 −> S6 0 0.5 0 0.5 1 S3 −> S6 0 0.5 0 0.5 1 S4 −> S6 0 0.5 0 0.5 1 S5 −> S6 Freq P D C GmAR−SS estimated PDC Significance threshold Figure 6.11: PDC of all the 6 estimated sources given by GmAR-SS (blue lines with circles) and the significance threshold determined using surrogate data (solid black line). 110 S1 S2 S3 S4 S5 S6 S2 S4 S3 S1 (a) Connectivity pattern of estimated sources (a) Connectivity pattern of simulated sources S1 S2 S3 S4 S5 S6 S2 S4 S3 S1 (a) Simulated sources S1 S2 S3 S4 S5 S6 S2 S4 S3 S1 (a) Connectivity pattern of estimated sources (a) Connectivity pattern of simulated sources S1 S2 S3 S4 S5 S6 S2 S4 S3 S1 (b) Estimated sources Figure 6.12: Connectivity patterns of simulated and GmAR-SS estimated sources get remained present (internally-guided). The two conditions were performed consecutively, each lasting for 15 s, and followed by a 15 s rest period. This sequence was performed twice at each of the three frequencies during one recording session. For PD subjects, the measure- ment was taken after a minimum of 12 hrs withdrawal of L-dopa medication. Calibration of the squeeze bulb to record maximum force was conducted at the beginning of each recording session. The EEG data were collected using an EEG cap (Quick-Cap, Compumedics, Texas, USA) with 19 electrodes using the International 10-20 system, referenced to linked mastoids. EEG data were sampled at 1000Hz using SynAmps2 amplifiers (NeuroScan, Compumedics, Texas, USA) and were later bandpass filtered off-line between 1 and 70Hz and down-sampled to 250Hz. Artifacts associated with eye blinks and muscular activity were removed using the Automated Artifact Removal in the EEGLAB Matlab Toolbox [82]. 6.4.2 Data Preprocessing For each subject, EEG signals recorded during the follow condition were concatenated. To reduce the computation load, the data are dimension reduced from 19 to 6 by dividing the electrodes into 6 regions and spatially averaging the EEG signals from the electrodes within each region [67]. These 6 regions were (see Fig. 6.14): Region 1 - (FP1, F7 and F3), Region 111 100 % 0 % % MVC 5 10 15 20 25 0 1 2 3 4 Time (sec) Ta rg et Figure 6.13: Squeezing task. The subject was instructed to follow the target bar (yellow) as close as possible. The force exerted by the subject was shown by the green bar. Oscillation of the target is shown in the bottom panel. 2 - (FP2, F4 and F8), Region 3 - (T7, C3, P7, P3), Region 4 - (T8, C4, P8, P4), Region 5 - (Fz, Cz, Pz) and Region 6 - (O1, O2). The spatial averaging is used here instead of PCA/PLS because we found in our preliminary analysis and prior study [109] that it helps minimize the effects of spurious activity measured from electrodes and yields more consistent connectivity results for this particular dataset. Nonetheless, it should be noted that in general, PCA-based dimension reduction techniques are more robust and should be used unless prior information is available. The proposed source separation algorithm was applied to the averaged signals of these 6 regions, yielding 6 source estimates. Let us denote the projection patterns from the sources to the 6 regions by Cregion and the projections from the sources to the 19 electrodes byCelec. To recoverCelec, linear regression was used: xi = s∗φ Ti +ε ∀ i= 1, . . . ,19 (6.24) where xi ∈RN×1 is a column vector representing the EEG signal from the i-th electrode and s∗ ∈ RN×6 is the time courses of all source estimates. The regression parameter, φ i ∈ R1×6, 112 FP1 FP2 Fz Cz Pz F7 F3 F4 F8 T7 P7 C3 P3 C4 T8 P4 P8 O1 O2 Region 1 Region 2 Region 3 Region 4 Region 5 Region 6 Figure 6.14: Six brain regions. The groupings of the electrodes are: Region 1 - (FP1, F7 and F3), Region 2 - (FP2, F4 and F8), Region 3 - (T7, C3, P7, P3), Region 4 - (T8, C4, P8, P4), Region 5 - (Fz, Cz, Pz) and Region 6 - (O1, O2). forms the i-th row of the mixing matrix estimate Celec. The resulting Celec is of dimension 19×6. The proposed source separation algorithm was applied to each subject separately with P= 6 as determined by AIC. The algorithm was run 10 times with random initial conditions for each subject and the solution with the highest likelihood was used as the final estimates. The clustering methods described in Section 6.2.3 was used to construct the global spatial filter. To illustrate the effects of experimental conditions on the source connectivity patterns for different subject groups, we concatenated the source signals from all subjects in each subject group (Normal/PD) for a given condition (follow/internally-guided) and fit a 6th-order mAR to the concatenated data. PDC estimates were computed over the frequency range between 1 and 40 Hz. The standard deviation of PDC estimates was determined using leave-one- subject-out technique. 6.4.3 EEG Results The clustering results including the dendrogram and the similarity matrix are shown in Fig. 6.15. Based on the tightness of the clusters, 20 clusters were identified. In particular, Cluster 8, 14, 18, 19 and 20 were identified as significant clusters with high inter-subject repeatability, and their scalp projection patterns, which are used to form the global spatial filter, are shown in 113 00.20.40.60.8 Cluster merging dissimilarity Dendrogram (linkage strategy: AL) Similarities between source projection patterns Red lines indicate estimate−clusters 18 417 10 19 9 14 15 1 8 11 13 7 3 20 16 212 6 5 Es tim at e− cl us te r l ab el 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 6.15: Dendrogram and similarity matrix. Fig. 6.17. The numbers of normal and PD subjects that contribute to a given cluster are as follows: Cluster 8: 4 out of 8 Normal, 3 out of 9 PD; Cluster 14: 7 out of 8 Normal, 9 out of 9 PD; Cluster 18: 4 out of 8 Normal, 8 out of 9 PD; Cluster 19: 3 out of 8 Normal, 3 out of 9 PD; Cluster 20: 7 out of 8 Normal, 6 out of 9 PD. The scalp projection patterns serve visually as a guide to interpreting what brain structures the isolated sources may represent. Cluster 8 was weighted mostly over the right sensorimo- tor cortex, while Cluster 14 was weighted over the left sensorimotor cortex. Cluster 19 was weighted over the left posterior parietal region, while Cluster 18 was weighted mostly occip- itally. Cluster 20 was widely distributed over the scalp with opposite weightings on either hemisphere. Cluster 20 may represent a deep source such as thalamic or basal ganglia struc- tures as have been previously described [165]. Before presenting the group analysis results, we provide an example of subject-level esti- mation results to illustrate the potential issue of inter-subject variability in real EEG analysis. The estimated scalp projection maps of a typical PD subject are shown in Fig. 6.16. Com- pared to the centrotypes of the significant clusters shown in Fig. 6.17, the scalp projection 114 (a) S1 (b) S2 (c) S3 (d) S4 (e) S5 (f) S6 Figure 6.16: Scalp projection maps of the estimated sources given by GmAR-SS from a typical PD subject. Blue dots represent positive weighting on the electrodes whereas red dots represent negative weighting. The radius of each circle is pro- portional to the voltage gain at the respective electrode. maps of the estimated sources S1, S2, S3 and S4 in Fig. 6.16 show a very high level of simi- larity to those of Cluster 8, 14, 18 and 20 in Fig. 6.17 which are the more dominant clusters in PD group. The estimated sources S5 and S6 cannot be related to any of the significant clusters, indicating the possible presence of subject-specific brain sources. PDC spectra of the source estimates from group spatial filtering are shown in Fig. 6.17. The majority of the PDC curves had peaks at 10Hz and/or 30Hz. In normal subjects, when they performed the internally-guided task, there was a distinct increase in 10Hz and 30Hz bands from Clusters 20 → 8 and 14 → 19. This was not seen with information flow in the opposite directions. In the follow condition, information flow was generally reduced, with the exception of a small increase in PDC in Clusters 18→ 20 at around 30Hz. In contrast, in PD subjects, the differences between the follow and internally-guided conditions were much less. Unlike controls, the increase in 10Hz and 30Hz bands from Clusters 20→ 8 was seen in the follow condition. There was also a large increase in 10Hz from Clusters 14 → 18 in 115 Cluster 8 10 20 30 40 0 0.2 0.4 Clust 14 −> Clust 8 10 20 30 40 0 0.2 0.4 Clust 18 −> Clust 8 10 20 30 40 0 0.2 0.4 Clust 19 −> Clust 8 10 20 30 40 0 0.2 0.4 Clust 20 −> Clust 8 10 20 30 40 0 0.2 0.4 Clust 8 −> Clust 14 Cluster 14 10 20 30 40 0 0.2 0.4 Clust 18 −> Clust 14 10 20 30 40 0 0.2 0.4 Clust 19 −> Clust 14 10 20 30 40 0 0.2 0.4 Clust 20 −> Clust 14 10 20 30 40 0 0.2 0.4 Clust 8 −> Clust 18 10 20 30 40 0 0.2 0.4 Clust 14 −> Clust 18 Cluster 18 10 20 30 40 0 0.2 0.4 Clust 19 −> Clust 18 10 20 30 40 0 0.2 0.4 Clust 20 −> Clust 18 10 20 30 40 0 0.2 0.4 Clust 8 −> Clust 19 10 20 30 40 0 0.2 0.4 Clust 14 −> Clust 19 10 20 30 40 0 0.2 0.4 Clust 18 −> Clust 19 Cluster 19 10 20 30 40 0 0.2 0.4 Clust 20 −> Clust 19 10 20 30 40 0 0.2 0.4 Clust 8 −> Clust 20 Frequency (Hz) P D C 10 20 30 40 0 0.2 0.4 Clust 14 −> Clust 20 10 20 30 40 0 0.2 0.4 Clust 18 −> Clust 20 10 20 30 40 0 0.2 0.4 Clust 19 −> Clust 20 Cluster 20 N−follow N−internal PD−follow PD−internal Figure 6.17: Diagonal subfigures: Scalp projection patterns corresponding to the centrotypes of Clusters 8, 14, 18, 19 and 20. Blue dots represent positive weighting on the electrodes whereas red dots represent negative weighting. The radius of each circle is proportional to the voltage gain at the respective electrode. Off-diagonal subfigures: PDC connectivity patterns between EEG sources derived from spatial filtering illustrating Normal vs. PD and Follow vs. Internally-guided. x-axis: frequency (Hz). y-axis: PDC value. 116 Cluster 8 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 14 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 18 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 19 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 20 Fr eq ue nc y (H z) Time (sec) 2 4 6 8 10 12 14 0 10 20 30 (a) Normal Subjects, Follow Condition Cluster 8 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 14 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 18 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 19 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 20 Fr eq ue nc y (H z) Time (sec) 2 4 6 8 10 12 14 0 10 20 30 (b) Normal Subjects, Internally-guided Condition Figure 6.18: Average spectrograms of the EEG sources of normal subjects derived from spatial filtering across all trials and all subjects of each subject group. Note the relative lack of apparent differences across tasks in groups compared to Fig. 6.17. 117 Cluster 8 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 14 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 18 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 19 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 20 Fr eq ue nc y (H z) Time (sec) 2 4 6 8 10 12 14 0 10 20 30 (a) PD Subjects, Follow Condition Cluster 8 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 14 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 18 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 19 Fr eq ue nc y (H z) 2 4 6 8 10 12 14 0 10 20 30 Cluster 20 Fr eq ue nc y (H z) Time (sec) 2 4 6 8 10 12 14 0 10 20 30 (b) PD Subjects, Internally- guided Condition Figure 6.19: Average spectrograms of the EEG sources of PD subjects. 118 the follow condition that was not seen in the 18 → 14 direction, suggesting an asymmetric information flow between Clusters 14 and 18. On the other hand, the frequency spectrograms of the extracted sources (see Fig. 6.18 and Fig. 6.19) did not show noticeable difference between groups and conditions. This result is consistent with the findings in our previous work [109] where we found that connectivity analysis was a more sensitive technique than spectral analysis in revealing the complex effects of motor tasks/disease. Our results suggest that PD subjects have a distinct inability to modulate activity from the sensorimotor cortex to other regions during internally-guided movements. This can be seen from examining the PDC from Cluster 14 to the other clusters (the second column in Fig. 6.17). In normal subjects, peaks at 10Hz (and to a lesser extent at 30Hz) can be seen from Cluster 14 during the internally-guided task. In contrast, PD subjects are unable to increase information flow during the internally-guided task and have an increase in activity from the left sensorimotor to the occipital clusters (14 → 18) during the follow condition. This inability to modulate activity from the sensorimotor cortex contralateral to the moving hand may be the result of “noisy” basal ganglia inputs to cortical structures [166]. The other striking difference between PD subjects and controls was the altered outflow of activity from Cluster 20 during internally-guided movements (the last column in Fig. 6.17). Again, normal subjects demonstrated peaks at 10Hz and 30Hz during internally-guided movements that were not replicated in PD subjects. The abnormalities that we detected in the PDC at 10Hz and 30Hz between many clusters in PD are compatible with prior studies using coherence [71], and examining STN recordings during DBS surgery. 6.5 Discussions One assumption of the proposed GmAR-SS method is that the innovation processes driv- ing the brain sources follow a non-Gaussian distribution in order to achieve good source separability (see Fig. 6.8). Although there is no physiologically motivated reasons for this choice, several previous studies have demonstrated meaningful non-Gaussian brain sources using ICA-based methods [167–169]. In our EEG analysis, we also found that identified 119 brain sources following super-Gaussian distributions, confirming the plausibility of the non- Gaussian assumption. As shown in Section 6.3, the proposed GmAR-SS yields better estimation performance compared to other ICA-based approaches. However, the performance improvement of GmAR- SS comes at the expense of increased computational cost which mainly is from two sources: one is the use of the gradient-ascend algorithm for finding solutions to Eqns. 6.13, 6.14 and 6.15, and the other comes from the use of the multi-start strategy for avoiding local maxima. In the case of ICA-based connectivity estimation techniques, such as MVAR-ICA, compu- tationally efficient gradient-descent ICA algorithms are available to allow fast estimation of model parameters. However, as noted in [170], the number of local optima in the ICA esti- mation can be affected by the distribution and signal-to-noise ratio of the data and thus the convergence to the global optimal is not always guaranteed in most real data application. While new ICA algorithms such as RADICAL [171] and MILCA [172] have been developed which use an exhaustive search over Jacobi angles to escape local optima, this comes at the expense of a substantial increases in the computation time for the estimation (see comparison in [173]). To circumvent these issues, gradient-descent based ICA algorithms, combined with the randomized multi-start technique, are often employed in order to obtain the global solu- tion. Given the computational efficiency of such ICA algorithms, the multi-start can be easily done within a short period of time, thus making MVAR-ICA a good choice when computation time is constrained. In the simulations, data dimension reduction was performed prior to model fitting and the correct dimension (J =M = 4) was used in all simulation runs since our evaluation schemes operate based on the assumption there is a one-to-one mapping between the simulated and estimated sources. However, the impact of dimensionality mismatch (i.e., J 6= M) on the estimation performance of existing connectivity estimation techniques has not yet been in- vestigated rigorously in the literature. This problem is particularly crucial for real data ap- plications since the true source dimensionality is generally unknown. A common practice used in the literature to determine source dimensionality is to choose a dimension such that 120 a certain percentage of variance (e.g., 99%) of the original data is retained after the dimen- sion reduction procedure [148, 158]. Clustering of the estimated sources based on their scalp projections maps or activation waveforms (see Section 6.2.3) may also be applied to help eliminate spurious or duplicate sources. Nonetheless, a more rigorous study on the effects of the dimensionality mismatch issue is warranted in the future work. 6.6 Conclusion To summarize, in this chapter we proposed a novel estimation method GmAR-SS that pro- vides robust estimations of the connectivity patterns between underlying brain sources using scalp EEG measurements. By incorporating a generalized mAR process into the previously proposed state-space mAR framework in [148], we are able to derive a maximum likelihood estimator that allows the mixing process and the AR model parameters to be jointly estimated from scalp EEG simultaneously. The proposed GmAR-SS technique, shown by simulations, yields considerably better estimation performance compared to conventional connectivity es- timation techniques under different simulation scenarios. Furthermore, the clustering tech- nique employed in this chapter resolves the source ambiguity problem commonly seen in applications of source separation techniques to group analysis. The real EEG results demon- strated that the proposed GmAR-SS approach is promising for discovering task-related brain source connectivity. 121 Chapter 7 A Multiblock PLS Model for Corticomuscular Interactions 4 7.1 Introduction One approach to assessing the final common pathway of the motor system is to simultane- ously assess cortical and muscle activity during sustained isometric muscle contraction. The most widely used method to study the interaction between the brain and the musculature or corticomuscular coupling is the magnitude-squared coherence (MSC). The MSC provides a normalized measure of correlation between two signals in the frequency domain. In mon- keys, beta range (13-50 Hz) coherence was found between electromyogram (EMG) and LFP recorded from the motor cortex during the hold phase of a precision grip task [174]. Simi- lar findings of beta-band corticomuscular coherence have also been demonstrated in humans using either the MEG [26, 27] or EEG [28] from the primary motor cortex during isometric contractions of weak to moderate strength. Alterations in cortico-muscular coherence have been demonstrated in motor-related diseases such as Parkinson’s disease [30]. 4A version of this chapter has been published as: J. Chiang, Z. J. Wang, and M. J. McKeown, “Multiblock PLS model for group corticomuscular activity analysis in Parkinson disease,” in Conference Record of the Forty Fourth Asilomar Conference on Signals, Systems and Computers (ASILOMAR’10), Nov. 2010, pp. 1375-1379. c© 2010 IEEE. http://dx.doi.org/10.1109/ACSSC.2010.5757759 122 Despite the wide use of MSC in studying corticomuscular coupling, it suffers from sev- eral limitations. The common formulation of coherence assumes pair-wise comparisons yet it has been suggested that the mapping between the brain and musculature during movement is many-to-many, as opposed to one-to-one [32]. Thus, coherence provides a limited view of the interactions between neural assemblies and muscles. In addition, EEG recordings typi- cally include electrical activity unrelated to motor performance, including eye and movement artifacts, and background brain activity. In fact, only a small fraction of the signal content is related to motor control [33]. Hence, when studying the relationship between EEG and EMG, the conventional approach of applying MSC directly to raw EEG and EMG signals typically yields a very low coherence value. Extensive statistical testing with e.g., bootstrap proce- dures, is required to demonstrate that the low-level coherence detectable is distinguishable from noise. Several multivariate data-driven approaches developed to analyze multi-modal data may circumvent some of the aforementioned concerns regarding MSC. In particular, techniques based on blind/semi-blind source separation have proven fruitful in multi-modal brain imag- ing analysis. Partial least square, which is widely applied in the chemical industry for process monitoring, exploits the covariation between predictor variables and response variables and finds a new set of latent components that maximally relate them [162]. The PLS method has been used to find spatial patterns of brain activity in fMRI data that are associated with behav- ioral measures [175] and to identify common temporal components between EEG and fMRI [176]. As such, the PLS method is well-suited for medical data analysis due to its ability to handle high dimensional and collinear data, where traditional regression-based techniques may break down unless extensive preprocessing is performed, such as variable selection and/or dimension reduction. Moreover, PLS is flexible enough to be extended to accom- modate multiple data sets in a multi-task/multi-subject scenario and also multiway data [177] such as EEG spectrogram which enables joint temporal, spectral and spatial decompositions as opposed to the traditional spatio-temporal decompositions. Although PLS has been used in various biomedical studies, the presence of inter-subject 123 variability is often overlooked or not dealt with explicitly when applying PLS to multi-subject data. One straightforward approach is to simply concatenate all subjects’ data along one of the data dimensions and perform the decomposition under the assumption that all subjects share common group patterns for the remaining data dimensions. This assumption is tenuous, as inter-subject variability may not be well modelled by simply pooling the data, especially in the presence of disease. Another approach is to apply PLS to each subject’s data separately, but then how to meaningfully match the derived components across subjects is unclear. To address this issue, we propose the utilization of a multiblock PLS (mbPLS) [178] approach, which is a two-level, multiple-block extension of the regular PLS. The mbPLS method groups data features of each subject into individual “subject blocks” at the sub-level, and aggregates the summary information from each subject block at the super-level, forming a “super-block” from which group-level components exhibiting high covariance between the predictor side and response side are extracted. In this hierarchical structure, subject-specific details are encapsulated at the sub-level and the group “consensus” at the super-level provides better interpretability for the multi-subject data. The mbPLS provides a trade-off between modeling simplicity and preservation of between-subject variation. In the current study, we apply the mbPLS method to concurrent EEG and EMG data collected from normal subjects and patients with PD performing a force tracking task. We exploit the ability of the proposed method to model 3-dimensional matrices and investigate time-frequency-connection relations using PDC to determine if ongoing connectivity between cortical regions corresponds with the ongoing EMG activity. We demonstrate that in both normal and PD subjects, the PDC between EEG channels is highly modulated in a manner matching the dynamic motor task performed and show significant correlations with EMG amplitude. 124 7.2 Methods 7.2.1 PLS Model and Estimation In this section, we give a brief overview of the basic two-block PLS and then introduce the mbPLS framework for group analysis. For notation, we use lower case letters, e.g., t , to denote column vectors and upper case letters, e.g., X , to denote two-dimensional matrices. Three-dimensional matrices are denoted by underlined upper case letters, e.g., X . Transposed vectors or matrices are indicated by a superscript T , e.g., tT . We use k (k = 1, . . . ,K) for indexing PLS components and b (b = 1, . . . ,B) for indexing subjects. The indices, t, f , c and m, are defined respectively with the following ranges: t = 1, . . . ,Nt, f = 1, . . . ,N f , c= 1, . . . ,Nc and m= 1, . . . ,Nm. Partial Least Square The two-block PLS model consists of two groups of variables, commonly preferred to as the predictor X and the response Y . The goal is to successively find orthogonal linear combina- tions of the predictor and response variables, known as predictor/response scores, that account for as much as possible of the covariation between X and Y . Specifically, PLS model can be represented by two outer relations that decompose the data blocks into sum of components X = TPT +E = K ∑ k=1 t kpTk +E , (7.1) Y = UQT +F = K ∑ k=1 ukqTk +F , (7.2) and an inner relation which ensures the maximal covariance between components uk = βkt k, k = 1, . . . ,K, (7.3) where K denotes the total number of PLS components. The vectors t k and uk are the scores of the k-th PLS component for X and Y , respectively, pk and qk are the associating normalized 125 SR PLS maximal covariance ResponsePredictor k1 r , k2 r , kB r , Rk k t T k1 p , T k2 p , T kB p , T sup w k u T sup q Sk k1 s , k2 s , kB s , T k1 q , T k2 q , T kB q , Y1 …Y2 YBX1 …X2 XB Figure 7.1: Multiblock PLS model. The mbPLS model consists of two levels: the super- level modeling the covariation between predictor and response variables and the sub-level collecting the details of individual blocks. loading vectors, and βk is the regression coefficient for the k-th component. The matrices E and F are the residual matrices. The number of components needed to describe the data blocks is typically determined based on the amount of variation remained in the residual matrices [162]. The estimation of the PLS model is done in a sequential fashion, component by compo- nent. The estimation starts with a random initialization of the response score u1. This vector is regressed on the predictor block X to give the block weight w1. The weight is normalized and multiplied with X to give the predictor score t1. For the response variables, the regres- sion is done similarly by regressing t1 on Y to yield a new u1. This is repeated until t1 and u1 converge to a predefined precision. After convergence, loading vector p1 is calculated by regressing t1 on X and data blocks are deflated by subtracting t1pT1 from X . The second pair of PLS components, orthogonal to the first, can be determined by repeating the iteration on the residual data blocks. The complete algorithm for estimating two-block PLS model is given in [162]. 126 Multiblock PLS The multiblock PLS extends the two-block PLS by dividing variables into multiple concep- tually meaningful data blocks, e.g., on a subject or trial basis. This blocking leads to two model levels: the super-level modeling the covariation between predictor and response vari- ables and the sub-level collecting features of individual blocks (Fig. 7.1). On the sub-level, the predictor variables and response variables of each subject/trial b are grouped into the pre- dictor block X b and the response blockY b, respectively. Each sub-level predictor block X b is iteratively decomposed into a sum of components. Each component k is an outer product of the predictor block score rb,k and the normalized loading vector pb,k. Similarly, the response block is decomposed into a set of response scores sb,k and loading vectors qb,k. The sub-level block scores are assembled to form the super blocks, Rk and Sk, on the super-level: Rk = [ r1,k,r2,k, . . . ,rB,k ] , (7.4) Sk = [ s1,k,s2,k, . . . ,sB,k ] . (7.5) These super blocks, representing a condensed subset of the sub-level variables, are treated as the predictor and response blocks as in the regular two-block PLS. From these, the predictor super score t k and response super score uk that best explain the variance in Rk and Sk while having maximal covariance between them are sought, along with the associating loading vectors, wsup and qsup. Specifically, the mbPLS model can be described as follows: 1. Sub-level decompositions: X b = K ∑ k=1 rb,kpTb,k+E b, (7.6) Y b = K ∑ k=1 sb,kqTb,k+F b, b= 1, . . . ,B. (7.7) 127 2. Super-level decompositions: Rk = [ r1,k,r2,k, . . . ,rB,k ] , (7.8) Rk = t kwTsup, (7.9) Sk = [ s1,k,s2,k, . . . ,sB,k ] , (7.10) Sk = ukqTsup. (7.11) 3. Inner relation between predictor and response blocks uk = βkt k. (7.12) The key advantage of mbPLS is that by organizing the data features into a hierarchical struc- ture, each block of subject features is allowed to have its own decomposition, thus allowing slight variations across subjects, but at the same time, the aggregate of the subject-level com- ponents are still constrained to achieve maximal covariance between the predictor blocks and response blocks. As in the regular two-block PLS, the estimation of mbPLS model is done sequentially, component by component [178]. The key difference in the estimation procedure is that in mbPLS, the regression is done at two levels: the predictor super score t k is regressed on the response super block Sk as well as on the sub-level blocksY b’s. This step allows the mbPLS to preserve the important features of regular two-block PLS while ensuring that maximal similarity between data blocks can be achieved. Beginning with a random initialization of the response super score u1, the sub-level predictor block scores, rb,1, and the corresponding loading vectors of the first PLS component are obtained by regressing u1 on X b individually for b= 1, . . . ,B. Subsequently, another regression is done on the super blocks R1 constituted by the new predictor block scores rb,1 to obtain the predictor super score t1. The response block scores and super scores are estimated similarly by regressing t1 on sub-level response data blocks and super block. This procedure is repeated until the super scores converge. The 128 algorithm of estimating mbPLS is detailed in Appendix A. Three-way Multiblock PLS In biomedical studies, some data are more meaningfully represented as a multiway array (or multidimensional matrix). For example, the EEG spectrogram can be expressed as a three- dimensional matrix of: (number of time points × number of frequency bins × number of channels). In this study, we extend the mbPLS model described in Section 7.2.1 to multiway arrays. The specific model considered in this work is a three-way mbPLS which consists of B three-dimensional predictor blocks of dimension (Nt×N f ×Nc) and B two-dimensional re- sponse blocks of dimension (Nt×Nm). Following the structure of a two-way mbPLS model, the three-way mbPLS model is expressed by two outer relations: Xb,t f c = K ∑ k=1 rb,k,tab,k, f pb,k,c+ eb,t f c, (7.13) Yb,tm = K ∑ k=1 sb,k,tqb,k,m+ fb,tm, (7.14) b= 1, . . . ,B, k = 1, . . . ,K, where Xb,t f c denotes the element (t, f ,c) in the three-dimensional matrix X b. The vectors rb,k = [rb,k,1,rb,k,2, . . . ,rb,k,Nt ]T and sb,k = [sb,k,1,sb,k,2, . . . ,sb,k,Nt ] are the score vectors, rep- resenting the temporal activations of the k-th component for X b and Y b, respectively. The vectors ab,k = [ab,k,1,ab,k,2, . . . ,ab,k,N f ]T and pb,k = [pb,k,1, pb,k,2, . . . , pb,k,Nc]T are the associ- ating normalized loading vectors ofX b which represent the spectral and spatial patterns of the k-th component, respectively. In addition, in order to obtain common spectral patterns across subjects for subsequent group analysis, we further assume that a1,k = a2,k = . . .= aB,k, k = 1, . . . ,K. (7.15) 129 As in the regular mbPLS, the relationship between the predictor blocks and response blocks is established by combining the sub-level block scores into super blocks Rk and Sk. From these, the super scores t k and uk that best explain the variance in Rk and Sk while having maximal covariance between them are sought. The estimation of three-way mbPLS model extends the three-way PLS algorithm proposed by [177] and is described in detail in Algorithm 1. Note that in Algorithm 1, the component index k is omitted from the variables for notational simplicity. 7.2.2 Application of mbPLS and EEG/EMG Feature Generation In this study, our goal is to find features in the EEG data whose temporal patterns are maxi- mally correlated with the EMG signals. In most existing studies, the corticomuscular analysis is performed directly on the raw EEG and EMG and typically yields very small correlation values. Nonetheless, with proper preprocessing steps, one can examine the EEG data for component(s) which have high correlation with the EMG signals. In this study, we exam- ined the relations between time-varying spectral features of the EEG signals, constituting the predictor variables of the mbPLS model, and amplitudes of the EMG signals and behavioral performance, constituting the response variables. A graphical representation of the EEG and EMG features and their decompositions are shown in Fig. 7.2. EEG Features Two types of EEG features were considered: 1) the EEG spectrogram which represents the time-varying spectral patterns of individual brain regions, and 2) the frequency-domain con- nectivity pattern between regions estimated using PDC. EEG Spectrogram Traditional approaches to studying EEG typically involve decomposing the EEG data into a set of frequency bands (alpha, beta, etc.) based on prior empirical knowl- edge. While averaging EEG activity over broad frequency bands can simplify subsequent 130 Algorithm 1 Calculation of the k-th Component of Three-way Multiblock PLS 1: Choose start u, e.g., u← the first principal component of some response blockY b. 2: repeat 3: Compute the auxiliary matrix Zb whose ( f ,c)-th element is given by zb, f c ← ∑Ntt=1Xb,t f cut , for f = 1, . . . ,N f ,c= 1, . . . ,Nc and b= 1, . . . ,B. 4: Concatenate the auxiliary matrices, Z ← [Z1Z2, . . . ,ZB]. The resulting Z is of dimen- sion (N f ×BNc) 5: Perform singular decomposition on Z : 6: [g,σ ,h]← SVD(Z ,1), where g and h are the first singular vectors of Z . 7: for b= 1 to B do 8: ab← g 9: pb← [ h(b−1)Nc+1,h(b−1)Nc+2, . . . ,h(b−1)Nc+Nc ]T 10: rb← [ rb,1, . . . ,rb,Nt ]T , where rb,t ← ∑N ff=1∑Ncc=1Xb,t f cab, f pb,c for t = 1, . . . ,Nt 11: end for 12: R← [r1,r2, . . . ,rB] 13: wsup←RTu/ ‖ uTu ‖, normalize wsup to ‖wsup ‖= 1 14: t ←Rwsup/wTsupwsup 15: for b= 1 to B do 16: qb←Y Tb t/ ‖ tTt ‖, normalize qb to ‖ qb ‖= 1 17: sb←Y bqb/qTbqb 18: end for 19: S← [s1,s2, . . . ,sB] 20: qsup← STt/ ‖ tTt ‖, normalize qsup to ‖ qsup ‖= 1 21: u← Sqsup/qTsupqsup 22: until t converges 23: Deflation by subtracting the effects of the first k components from each data block, X b andY b for b= 1, . . . ,B: 24: for b= 1 to B do 25: Gb← (T T1:kT 1:k)−1T T1:kub, where T 1:k = [t1, . . . ,t k]. 26: Xb,t f c← Xb,t f c− rb,tab, f pb,c for t = 1, . . . ,Nt, f = 1, . . . ,N f and c= 1, . . . ,Nc 27: Y b←Y b−T 1:kGbqTb 28: end for 29: k← k+ 1. The next pair of PLS components is found by repeating from Step 1 on the deflated data blocks. 131 EEG Features XbTime Frequency EEG channel = ∑ k=1 K Normal: Spatial Pattern PD: Spatial Pattern kbp , kba , kbr , (a) EEG Spectrograms 10 20 30 40 50Hz XbTime Frequency Connectivity between EEG channels = kbp , kba , kbr , (b) EEG Frequency-domain Connectivity 5 1 0 1 5 2 0 2 5 Se c 5 1 0 1 5 2 0 2 5 Se c 10 20 30 40 50 Hz ∑ k=1 K EMG Features Time [ EMG channels, Force ] = kbq , kbs , EMG Amplitude and Force Output Yb 5 1 0 1 5 2 0 2 5 Se c ∑ k=1 K Figure 7.2: Top - EEG features: (a) EEG spectrogram and its decomposition. (b) EEG frequency-domain connectivity and its decom- position. An EEG data block is decomposed into K components where each component k is an outer product of temporal (rb,k), spectral (ab,k) and spatial (pb,k) patterns. Bottom - EMG features: EMG data block consisting of EMG amplitude and force output, and its decomposition. An EMG data block is decomposed into K components where each component k is an outer product of temporal pattern (sb,k) and loading vector (qb,k). 132 analysis, it inevitably sacrifices the frequency resolution. In this study, instead of using band- limited EEG waveforms as an EEG feature for PLS, we propose using the time-varying EEG spectrogram which provides a simultaneous time-frequency-space decomposition. The EEG spectrogram was computed using the short-time Fourier transform (‘spectro- gram’ function in Matlab) for each brain region over the frequencies 3-50 Hz with a Ham- ming window of length 300 time points and a 95% overlap. The EEG spectrograms of each subject were organized as a three-dimensional matrix (Nt×N f ×Nc) with rows being time (t), columns being frequency ( f ) and the third dimension being brain region (c). The EEG spec- trograms were then decomposed by mbPLS into a sum of components, each with a specific temporal, spectral and spatial pattern. The spectral and spatial patterns allow the identifica- tion of brain oscillatory rhythms and their distributions over the scalp that give a temporal pattern having maximal covariance with the corresponding EMG temporal pattern. Frequency-Domain Connectivity Pattern: Partial Directed Coherence In addition to estab- lishing the relation between individual brain regions and active muscles, we were also inter- ested in identifying the frequency-dependent cortical couplings that were maximally mod- ulated by the motor task. To provide a spectral quantification of coupling between EEG signals, we used PDC [13]. PDC has been applied in various EEG studies as it provides the ability to distinguish between direct and indirect causal relations between signals and it has shown promising results in assessing the EEG changes in Parkinson’s disease [109]. The calculation of PDC starts with fitting a P-th order mAR model to the observed EEG signals. The mAR model parameters are estimated using the ARfit toolbox in Matlab [164]. Once the mAR coefficients are determined, they are transformed into the frequency domain via the Fourier transform, and the PDC between each pair of nodes (or regions of nodes) is calculated according to Eqn. 2.7 (refer to Section 2.2.2 for a detailed description). The mAR model assumes stationarity of the observed signals. In order to investigate the temporal evolution of connectivity patterns during dynamic movements, we used a moving window of length 300 time points with 95% overlap and computed PDC independently for 133 each time window over the frequencies 3-50 Hz. The mAR model order was determined by Schwarz’s Bayesian Criterion, and in this study, a 5-th order (P = 5) mAR model was used. The PDC values of each subject were organized as a three-dimensional matrix with rows being time (t), columns being frequency ( f ) and the third dimension being directional connection between brain regions (c). EMG Features An EMG signal can be considered as a zero-mean, band-limited and wide-sense stationary stochastic process modulated by the EMG amplitude [179]. The surface EMG amplitude is a reflection of overall muscle activity of individual underlying muscle fibers [179]. Although different techniques have been proposed for accurate amplitude estimation, we used the sim- ple root-mean-square (RMS) method, with a moving window of length 300 time points and 95% overlap applied. Note that the window size used here is same as that in EEG spectro- gram calculation, to ensure that the EEG and EMG features were temporally aligned. Sim- ilarly, the force output from the subject was resampled to ensure its proper alignment with EMG amplitude. The response variables, EMG amplitude and force output, were organized as a two-dimensional matrix (Nt ×Nm) with rows being time (t) and columns being EMG channels and force output (m). 7.2.3 Assessment of Significance To determine the statistical significance of the correlation between the temporal patterns of EEG and EMG, we used the non-parametric permutation test [180] in which the temporal correlation between EEG and EMG was removed by permuting the temporal order of the EEG featuresX b uniformly for all subjects while leaving the EMG featuresY b intact. For this study, 100 random permutations were generated and mbPLS was applied to each of these per- mutations. The correlation coefficients between the permuted EEG features and unchanged EMG features were then computed to form the empirical null distribution. The p-value of the original EEG/EMG correlation can be calculated from the null distribution as the propor- 134 tion of sampled permutations whose correlations are greater than or equal to the correlation between the original EEG and EMG features. In addition, we also want to determine which EEG features in the spatial patterns have significant contributions to the EEG/EMG temporal correlation. This is done by identifying EEG features that have weights statistically significant different from zero. In particular, we focus on the PLS components that have significant temporal EEG/EMG correlations. The significance level of the EEG spatial patterns was determined using the jackknife resampling method [180] from which delete-one jackknife samples are generated by leaving out time points from both X b and Y b, for b = 1, . . . ,B, one at a time and applying mbPLS to the remaining time points. Let p̂b( j) denote the EEG spatial pattern of data block b estimated from the delete-one jackknife sample that has time point j removed. The j-th jackknife pseudo pattern for the k-th component is given as p̂∗b,k( j) = Ntpb,k− (Nt−1)p̂b,k( j), (7.16) where Nt denotes the total number of time points in X b and pb is the spatial pattern of the full data. The mean, p̂∗b,k, and standard deviation, SE(p̂ ∗ b,k), of the jackknife pseudo patterns are calculated as p̂∗b,k = 1 Nt Nt ∑ j=1 p̂∗b,k( j), (7.17) SE(p̂∗b,k) = √√√√ 1 Nt Nt ∑ j=1 (p̂b,k( j)− p̂∗b,k)2. (7.18) From these, t-statistic for each element of pb,k is computed and forms the corresponding t- spatial pattern for pb,k. To allow for easier visualization and comparison of spatial patterns between different subject groups (PD subjects vs. normal subjects), for each significant PLS component, we concatenate the t-spatial patterns of subjects within a subject group horizon- tally and perform PCA on the concatenated t-spatial patterns. The fist principal component is 135 used to represent the group-level spatial pattern. 7.3 Experimental Data Data used in the study are the same as those used in Chapter 6. Detailed description of subjects and experimental procedure are presented in Section 6.4.1. In this study, we focus on the time periods when subjects were squeezing at a frequency of 0.4 Hz under the follow condition - both the visual target and visual feedback of the force generated were shown. We next describe the acquisition and preprocessing of EEG and EMG data. 7.3.1 EEG Data Acquisition and Preprocessing The EEG data were collected using an EEG cap (Quick-Cap, Compumedics, Texas, USA) with 19 electrodes using the International 10-20 system, referenced to linked mastoids. There EEG data were sampled at 1000 Hz using SynAmps2 amplifiers (NeuroScan, Compumedics, Texas, USA). A surface electrode on the tip of the nose was used as ground. Ocular movement artifacts were measured using surface electrodes placed above and below the eyes (Xltek, Ontario, Canada). Data were later bandpass filtered off-line between 1 and 70Hz and down- sampled to 250 Hz. Artifacts associated with eye blinks and muscular activity were removed using the Automated Artifact Removal in the EEGLAB Matlab Toolbox [82]. To simplify analysis, the EEG signals were divided into 5 regions: Fronto-Central (FCentral) - FP1, FP2, F7, F3, Fz, F4 and F8, Left Sensorimotor (LSM) - T7, C3, P7 and P3, Right Sensorimotor (RSM) - T8, C4, P8 and P4), Central - Cz and Pz and Occipital - O1 and O2 (see Fig. 3.2 for groupings of electrodes). The raw time courses of the electrodes within each region were averaged as the overall activity of the region, and the averaged time course were then zero-meaned and normalized to unit variance. For subsequent analysis, data collected during the squeezing periods con- catenated in time into a single matrix for each individual subject. Data from the rest periods are excluded from the analysis. 136 7.3.2 EMG Data Acquisition and Preprocessing The EMG signal of the wrist flexor compartment was recorded using self-adhesive, silver, silver-chloride pellet surface electrodes with 7 mm diameter. A bipolar montage was used with a fixed inter-electrode distance of 30 mm. The surface EMG signal was simultaneously collected together with the EEG signals and was amplified and sampled at 1000 Hz. To be consistent with the EEG preprocessing, the EMG signal was down-sampled off-line to 250 Hz and only the squeezing periods are used for subsequent analysis. 7.4 Results An example of the coherence spectra between the raw EEG signals of the 5 brain regions and the raw EMG signal of the right wrist flexor of a typical normal subject based on a non- overlapping Hamming window of length 300 points (1.2s) is shown in Fig. 7.3. The 95% confidence limit for the zero coherence is indicated by the dashed line. The 95% confidence limit is determined by [181] CL(α = 0.95) = 1− (1−α) 1L−1 , (7.19) where L is the number of windowed segments. Estimated coherence values above the CL were regarded as significant. Notice the small values of the estimated coherence, with peaks at ∼10Hz and ∼30Hz. The spectrograms of raw EEG signals of the same normal subject are shown in Fig. 7.4. None of the 5 regions shows clear modulation of activity by the motor task. The proposed mbPLS method was applied to concurrent EEG and EMG data of all 10 normal and 9 PD subjects simultaneously, which gave 19 predictor blocks and 19 response blocks in total. The joint modeling of normal and PD data allowed the identification of common temporal and spectral patterns. The spatial patterns are free to vary across subjects, from which one can identify brain regions that are differentially recruited by PD during a motor task. 137 5 10 15 20 25 30 35 40 45 50 0 0.2 0.4 0.6 0.8 1 Co he re nc e Frequency (Hz) Figure 7.3: Coherence between the raw EEG signal of the 5 region and the EMG signal of a typical normal subject. The dashed line represents the 95% confidence limit. Figure 7.4: Spectrograms of raw EEG signals. 138 5 10 15 20 25 30 35 40 45 50 0 0.1 0.2 0.3 0.4 Spectral Signature Frequency (Hz) Component 1, CC = 0.45 (p = 0.002) 5 10 15 20 25 −4 −2 0 2 4 Time (sec) Temporal Signature EEG EMG Target Normal: Spatial Pattern PD: Spatial Pattern Figure 7.5: Component 1 of PLS decomposition with EEG Spectrograms as predictors and EMG amplitude as responses. Top panel: Temporal patterns of EEG (solid) and EMG (dashed) super scores. The oscillation of the target bar is also shown (grey, solid). Middle panel: Spectral pattern. Bottom panel: Spatial patterns of normal subjects (left) and PD subjects (right). The radius of each circle is made proportional to the weighting of the respective electrode in the group EEG spatial pattern. 139 5 10 15 20 25 30 35 40 45 50 0 0.1 0.2 0.3 0.4 Spectral Signature Frequency (Hz) Component 7, CC = 0.21 (p = 0.114) 5 10 15 20 25 −4 −2 0 2 4 Time (sec) Temporal Signature EEG EMG Target Normal: Spatial Pattern PD: Spatial Pattern Figure 7.6: Component 7 of PLS decomposition with EEG Spectrograms as predictors and EMG amplitude as responses. 140 7.4.1 EEG Spectrogram Components Corresponding to EMG and Behavioral Data Using the permutation test, one significant component (p< 0.05) from the EEG spectrogram was identified. The temporal and spectral patterns for Component 1 are shown in Fig. 7.5. The correlation coefficient between the EEG t1 and EMG u1 super scores was 0.45, corresponding to p = 0.002. The spectral patterns shows peaks around 5-15 Hz which were associated with the theta and alpha rhythms. The group-level EEG spatial activation patterns of normal and PD subjects given by PCA are shown in the bottom panel of Fig. 7.5. Normal and PD subjects demonstrated similar spatial distributions in which higher values are located in the left sensorimotor area, contralateral to the performing hand. The temporal and spectral patterns for a component with a clear peak at around 30Hz (Component 7) are shown in Fig. 7.6. The correlation coefficient between EEG and EMG su- per scores is 0.21, which was not significant (p= 0.114). In normal subjects, this components is associated with activity in frontal and right sensorimotor regions, while in PD subjects, the activity is seen more in the frontocentral region. 7.4.2 PDC Components Corresponding to EMG and Behavioral Data Using the permutation test, three significant components from decomposition of the PDC (p < 0.05) were identified. The temporal and spatial patterns for Component 1 are shown in Fig. 7.7. The correlation coefficient between the EEG t1 and EMG u1 super scores was 0.58 (p= 0.0084). The associated spectral pattern demonstrated higher power in the beta and gamma bands. The group-level EEG spatial activation patterns of normal and PD subjects given by PCA are shown in the bottom panel of Fig. 7.7. For easier interpretation, only con- nections whose weights were greater than 80% of the maximum value are shown. In normal subjects, this component was dominated by the connection Central→ FCentral, whereas in PD subjects, it was dominated by connections Occipital→ FCentral and RSM→ FCentral. The temporal and spatial patterns for Component 2 are shown in Fig. 7.8. The correlation coefficient between EEG and EMG super scores is 0.54 (p= 0.0179) and had a broad spectral 141 pattern, with higher power concentrated in the lower beta range (15-25 Hz). In both normal and PD subjects, this component was dominated by connections going into the Occipital and RSM regions The temporal and spatial patterns for another significant component shown in Fig. 7.9. The correlation coefficient between EEG and EMG super scores is 0.43 (p = 0.0337). The spectral pattern of this component was concentrated between 10 and 25 Hz. In normal sub- jects, this component was based largely on connections RSM→ LSM and Central→ Occip- ital, whereas in PD subjects, it was based on connection LSM→ FCentral. 7.5 Discussion We have demonstrated that the proposed mbPLS method provides a framework for multi- modal, multi-subject analysis. It ensures similar decompositions across subjects by constrain- ing the block scores of every subject to have maximal covariance with the common group- level super score, yet details pertaining to individual subjects are preserved. When applied to EEG/EMG data during performance of the dynamic motor task of sinusoidal squeezing, we found results that were not apparent in any of the raw spectrograms (Fig. 7.4). Use of MSC detected a small but significant peak at ∼10Hz and ∼30Hz (Fig. 7.3). A component with a similar peak at ∼30Hz approached statistical significance (Fig. 7.6), but was concentrated more over frontal regions. The most interesting part of our results was that the PDC measures (Figs. 7.7-7.9) cor- related much better with ongoing EMG activity than activity at a single region (Figs. 7.5 and 7.6). At first it may seem counter-intuitive that the interaction between cortical regions would correlate more closely with EMG, since anatomically considerations would argue that the motor cortex would demonstrate the greatest correspondence. However, as previously noted, only a small fraction of the EEG signal at a given region is directly related to motor control [33]. Examining the interactions between regions may therefore constrain looking at the fraction of the EEG important for performance of this dynamic, visually-related task. We found that connections to and from occipital regions were prominent in PD subjects 142 5 10 15 20 25 30 35 40 45 50 0.08 0.1 0.12 Spectral Signature Frequency (Hz) Component 1, CC = 0.58 (p = 0.0084) 0 5 10 15 20 25 −4 −2 0 2 4 Time (sec) Temporal Signature EEG EMG Target Normal: Spatial Pattern PD: Spatial Pattern Figure 7.7: Component 1 of mbPLS decomposition with EEG PDC as predictors and EMG amplitude as responses. Top panel: Temporal patterns of EEG (solid) and EMG (dashed) super scores. The oscillation of the target bar is also shown (grey, solid). Middle panel: Spectral pattern. Bottom panel: Spatial patterns of normal subjects (left) and PD subjects (right). The width of each arrow reflects the relative weighting of the respective connections in the group EEG spatial pattern. 143 5 10 15 20 25 30 35 40 45 50 0.08 0.1 0.12 Spectral Signature Frequency (Hz) Component 2, CC = 0.54 (p = 0.0179) 0 5 10 15 20 25 −4 −2 0 2 4 Time (sec) Temporal Signature EEG EMG Target Normal: Spatial Pattern PD: Spatial Pattern Figure 7.8: Component 2 of mbPLS decomposition with EEG PDC as predictors and EMG amplitude as responses. 144 5 10 15 20 25 30 35 40 45 50 0 0.05 0.1 0.15 0.2 Spectral Signature Frequency (Hz) Component 8, CC = 0.43 (p = 0.0337) 0 5 10 15 20 25 −4 −2 0 2 4 Time (sec) Temporal Signature EEG EMG Target Normal: Spatial Pattern PD: Spatial Pattern Figure 7.9: Component 8 of mbPLS decomposition with EEG PDC as predictors and EMG amplitude as responses. 145 (Fig. 7.7 and 7.8). Relative to controls, PD subjects heavily on visual cues for the initi- ation [182] and ongoing control of movement [183]. In addition, the increased intra and inter-hemispheric connectivity we observed in the PD subjects (Fig. 7.8) is consistent that previously described with MEG during the resting state [4]. While several data driven techniques have been proposed for multi-modal data analysis, they still have limitations that may limit their suitability for medical applications. Techniques related to ICA and PCA relies on independence or orthogonality assumption which is often too stringent and has no obvious physiological justification. Most of the proposed techniques are inadequate for handling multi-subject data sets unless the data are combined a priori by grand-averaging or concatenation across subjects. However, both procedures rely on the assumption that variability between subject is negligible which, however, is rarely true for real data, especially in disease states such as Parkinson’s disease we examined here. In this study, we proposed the multiblock PLS framework to address both of these issues. The mbPLS provides a unique decomposition of data without assuming independence between components. Furthermore, with the hierarchical two-level design, mbPLS is able to find a middle ground between grand-averaging across subjects and doing decomposition on a subject-by-subject basis. At the sub-level, data blocking allows the preservation of subject specificity while at the super-level, aggregation of sub-level results allows the identification of common features across subjects. We note that while the proposed technique was applied to joint EEG and EMG coupling analysis, it is sufficiently general that it can be also applied to other modalities of neurophys- iological signals including MEG, fMRI and kinematic data. 7.6 Conclusion EEG and simultaneously-recorded EMG data are a means to assess integrity of the functional connection between the cortex and the muscle during movement. EMG-EEG coupling is typically assessed with pair-wise squared coherence, resulting in a small, but statistically- significant coherence between a single EEG and a single EMG channel. However, a means 146 to combine results across subjects is not straightforward with this approach because the exact frequency of maximal EEG-EMG coupling may vary between individuals, and it emphasizes the role of an individual locus in the brain in driving the muscle activity, when interactions between brain regions may in fact be more influential on ongoing EMG activity. To deal with these issues, we implemented a mbPLS procedure, previously proposed in chemical applications, which incorporates a hierarchical structure into the ordinary two-block PLS often used in neuroimaging studies. In the current implementation, each subject’s data fea- tures are collected in individual data blocks on a sub-level, while simultaneously aggregating the sub-level information to obtain a super-level group “consensus”. We further extended the mbPLS model to include 3-dimensional matrices: time-frequency-EEG channel and a time-frequency-connection utilizing PDC. We applied the proposed method to concurrent EEG and EMG data collected from ten normal subjects and nine patients with mild-moderate Parkinson’s disease performing a dynamic motor task - that of sinusoidal squeezing. The re- sults demonstrate that connections between EEG electrodes, rather than activity at individual electrodes, correspond more closely to ongoing EMG activity. In PD subjects, there was en- hanced connectivity to and from occipital regions, likely related to the previously-described enhanced use of visual information during motor performance in this group. The proposed mbPLS framework is a promising technique for performing multi-subject, multi-modal data analysis and it allows for robust group inferences even in the face of large inter-subject vari- ability. 147 Chapter 8 Conclusion This chapter will summarize and discuss the main findings of the thesis and provide sugges- tions for future research directions. 8.1 Contributions In this thesis, we developed new signal processing methods for inferring brain connectiv- ity from EEG recordings. The proposed connectivity analysis techniques provide enhance- ments to the previously proposed mAR model-based connectivity measures to address some of the key challenges present in real EEG studies such as volume conduction, low task-related signal-to-noise ratio, statistical instability in model estimation and group inference for multi- subject data. These techniques were applied to EEG data collected from healthy subjects and PD patients, including ‘on’ medication and ‘off’ medication states during the performance of motor tasks to study effects of disease and medication on connectivity patterns. The main contributions and findings of this thesis are summarized as follows. In Chapter 2, we investigated the statistical instability associated with estimating param- eters of a full mAR model from high-dimensional EEG signals and proposed a sparse-mAR based connectivity measure, called sparse PDC, where the sparsity is achieved by imposing a LASSO penalty on the coefficient matrices in mAR estimation. We show by simulation that the use of LASSO penalty facilitates automatic variable selection and yielded a more 148 stable and accurate estimate of connectivity compared to the traditional, non-sparse approach when data samples are limited. As the sparse formulation of PDC is more closely matching the sparse nature of brain connectivity and naturally addresses the concern of model order selection, it provides a useful tool for EEG connectivity analysis. In Chapter 3, we applied the sparse PDC technique to EEG data of normal and PD sub- jects, including both before and after Levodopa administration, while performing a visu- ally guided joystick task. The sparse-PDC analysis revealed multifaceted, regionally and directionally-dependent alterations of connectivity in PD subjects withdrawn from medica- tion during both movement preparation and execution. Connectivity was particularly altered posteriorly, suggesting abnormalities in visual and visuo-motor processing in PD. For com- parison, we utilized traditional Fourier analysis to evaluate task-dependent frequency spectra modulation in these same regions. While there were some overall spectral differences be- tween subject groups, the differences between regions tended to be relatively modest. In contrast, directional connectivity analysis demonstrated distinctive regional difference within groups. In particular, connectivity measures in the alpha, beta and low gamma frequency ranges showed high correlation with motor UPDRS in PD subjects withdrawn from medi- cation, suggesting the possibility of using connectivity as a potential biomarker for patho- logical brain oscillations in PD. When comparing connectivity patterns between on- and off- medication states of PD subjects, we observed that Levodopa medication only partially re- stored connectivity, and in some cases resulted in further exacerbation of abnormalities. Our results suggest that PD is associated with significant alterations in connectivity between brain regions and that these changes can be non-invasively detected in the EEG using the PDC method. Motivated by the observation of the association between connectivity and clinical mea- sure of motor function in PD subjects in Chapter 3, we proposed a LASSO method in Chapter 4 to select PDC-derived connection/frequency band combinations that best predict UPDRS scores. The predictive performance of the proposed method was assessed using the leave- one-subject-out cross validation procedure. Our results showed that particular combinations 149 of directional connectivity strengths were able to clearly differentiate unmedicated PD sub- jects from normal subjects and accurately predict clinical UPDRS scores of PD subjects. The frequencies selected by the LASSO were not restricted to the beta band as typically seen in PD literature but included both alpha and beta. The predicted UPDRS values in medicated PD subjects were normalized to a similar level as normal subjects. These observations cor- roborate the plausibility of using EEG-derived connectivity as a biomarker for evaluating PD deficits. In Chapters 3 and 4, inference of group-level connectivity was done by applying sparse PDC to each subject separately with no prior assumptions on the structural patterns across subjects and the results were later integrated across subjects using statistical methods, such averaging or ANOVA. In Chapter 5, we extended the sparse mAR model developed in Chapter 2 by imposing a common group structure across subjects. This is achieved using a data-driven approach by incorporating the Group LASSO method into the joint estimation of mAR coef- ficients of all subjects. The gLASSO enforces the mAR models of the subjects within a group to be structurally identical across subjects, but the connection coefficients are allowed to vary from subject to subject. The proposed group sparse mAR model represents a compromise be- tween the two extremes of conventional group analysis approaches as it allows identification of common group structures while preserving the subject-specificity. To address the issue of volume conduction, in Chapter 6, we proposed a state-space gener- alized mAR framework that jointly models the instantaneous mixing effects of volume con- duction and the causal relationships between underlying brain sources. Different from the existing two-stage source estimation approaches, we developed a maximum likelihood esti- mator allowing the the mixing process and mAR model parameters to be jointly estimated si- multaneously. The proposed GmAR source separation technique was verified using simulated EEG data based on a single-shell spherical head model with varying levels of signal-to-noise ratio, signal-to-background noise ratio, data lengths and source Gaussianity. The proposed technique demonstrated improved estimation accuracies compared to conventional two-stage connectivity estimation approaches. The GmAR-SS technique was applied together with 150 the clustering technique to real EEG data of normal and PD subjects performing a force- tracking task. The results demonstrated that the proposed group analysis framework resolves the source ambiguity issue commonly seen in group analysis and is effective in detecting differences between subject groups and experimental conditions. In Chapters 2-6, we focused on establishing the causal relationships between brain re- gions during motor performance. To assess the final pathway of the motor system, in Chapter 7, we expanded our connectivity analysis to the cortico-muscular level to look at interac- tions between the brain and musculature. Specifically, we proposed modeling simultane- ously recorded EEG and EMG signals using a multiblock PLS framework. The multiblock PLS framework provides a natural way to model multi-subject, multi-modal data. The PLS method exploits the covariation between two different sets of data variables (one being EEG and the other being EMG) and finds a new set of latent components that maximally relate them. The multiblock construction of PLS allows data to be divided into separate data blocks on a per-subject basis where each block can have its own decomposition, while at the same time, the aggregates of the subject-level EEG and EMG components are still constrained to achieve maximal covariance, allowing identification of common features across subjects. We further extended the mbPLS model to include 3-dimensional matrices: time-frequency-EEG channel and time-frequency-connection based on PDC. We applied the proposed method to concurrent EEG and EMG data collected from normal PD patients performing a sinusoidal squeezing task. The results showed that connections between EEG electrodes, rather than activity at individual electrodes, correspond more closely to ongoing EMG activity. 8.2 Future Research Directions Studying brain connectivity in the source domain provides a plausible solution to the volume conduction effect inherent in EEG analysis. However, the presence of inter-subject variabil- ity, which potentially leads to different brain sources being recruited and variations in the spatial distribution and strength of source activation across subjects, has made group infer- ence in source separation a challenging issue and prevented it from being widely applied in 151 real EEG studies. To address this, the approach we took in Chapter 6 was to learn the source separately for each subject and then perform group inference on the subject-specific sources. This was done by clustering the extracted sources based on their spatial distributions and se- lecting the ones that repeatedly appeared across subjects as the group representatives. In con- trast to the conventional approach of applying source separation to averaged or concatenated multi-subject data which in effect ignores the inter-subject variability, the clustering method represents a more flexible approach for dealing with inter-subject variability. However, the results can be sensitive to the choice of similarity measures and the number of clusters used to represents the entire set of sources. A potentially more robust method to perform group analysis in source separation is to incorporate a hierarchical model into the current GmAR-SS framework. Using a similar concept as in parametric empirical Bayesian framework [184], we assume that individual subjects’ connectivity patterns are generated probabilistically from a common, group structure. As a simple example, for each connection in the subject-level structure, its connectivity can be represented by a Bernoulli random variable whose parameter p is specified in the group structure based on anatomically motivated prior knowledge. When p= 1, the connection must be present in every subject and when p= 0, the connection can- not exist in any subject. Furthermore, the strengths of connections between sources (i.e., the mAR parameters) are allowed to vary across subjects. This hierarchical model can be further extended to include more hyper-parameters that may correspond to different experimental factors, such as disease state, before/after medication, etc. In addition, motivated by our findings in Chapter 7, another possible extension to the GmAR-SS framework is to include task-related outputs, such as EMG signals for motor tasks, or prior anatomical knowledge to constrain temporal or spatial patterns of the source compo- nents, allowing for more accurate identification of task-related sources. This in effect helps minimizing variability across subjects and simplifies subsequent group inference and inter- pretation of the results. Another possible research direction is to incorporate different decomposition techniques into the mbPLS framework. The original formulation of PLS utilizes PCA for decomposi- 152 tion. PCA decomposes a two-dimensional matrix into orthogonal subspaces. PCA is compu- tationally efficient but the orthogonality assumption is less intuitive from a neurobiological perspective. A possible alternative is to incorporate ICA into the PLS/mbPLS framework for data decomposition. The key difference is that in the PCA-based PLS, each successively extracted latent variable (LV) is constrained to be orthogonal to the preceding LVs whereas in ICA-based PLS, each LV is constrained to be statistically independent from the preced- ing LVs. In either case, the predictor LVs and response LVs are optimized to be maximally correlated. The ICA-based PLS is related to the concept of semi-blind ICA or constrained ICA which constrains the time courses of the estimated ICs to be ‘close’ to the experiment paradigm time course. They differ in the sense that in semi-blind ICA [185], only one data matrix (e.g., EEG) is being decomposed whereas in ICA-based PLS, both the predictor and responses variables (e.g., EEG and EMG) are decomposed simultaneously. The PLS/mbPLS framework can also be extended to model more than two sets (one being predictor and the other being response) of data variables. This can be useful in analyzing multiple simultaneously recorded physiological signals such as fMRI-EEG-EMG data fusion. Different approaches to enforce maximal joint correlations among multiple data variables, compared to maximizing pair-wise correlation between every possible pair of variables, in the PLS/mbPLS framework can be explored. 153 Bibliography [1] M. Belmonte, G. Allen, A. Beckel-Mitchener, L. Boulanger, R. Carper, and S. Webb, “Autism and abnormal development of brain connectivity,” The Journal of Neuroscience, vol. 24, no. 42, pp. 9228–9231, 2004. [2] A. Leuchter, T. Newton, I. Cook, D. Walter, S. Rosenberg-Thompson, and P. Lachenbruch, “Changes in brain functional connectivity in Alzheimer-type and multi-infarct dementia,” Brain, vol. 115, no. 5, pp. 1543–1561, 1992. [3] S. Farmer, “Neural rhythms in Parkinson’s disease,” Brain, vol. 125, no. 6, pp. 1175–1176, 2002. [4] D. Stoffers, J. Bosboom, J. Deijen, E. Wolters, C. Stam, and H. Berendse, “Increased cortico-cortical functional connectivity in early-stage Parkinson’s disease: An MEG study,” NeuroImage, vol. 41, no. 2, pp. 212–222, 2008. [5] K. Friston, “Functional and effective connectivity in neuroimaging: a synthesis,” Human Brain Mapping, vol. 2, no. 1-2, pp. 56–78, 1994. [6] B. Horwitz, “The elusive concept of brain connectivity,” NeuroImage, vol. 19, no. 2, pp. 466–470, 2003. [7] J. Schoffelen and J. Gross, “Source connectivity analysis with MEG and EEG,” Human Brain Mapping, vol. 30, no. 6, pp. 1857–1865, 2009. [8] S. Makeig, A. Bell, T. Jung, and T. Sejnowski, “Independent component analysis of electroencephalographic data,” in Advances in Neural Information Processing Systems, 1996, pp. 145–151. [9] E. Pereda, R. Quiroga, and J. Bhattacharya, “Nonlinear multivariate analysis of neurophysiological signals,” Progress in Neurobiology, vol. 77, no. 1-2, pp. 1–37, 2005. [10] V. Sakkalis, “Review of advanced techniques for the estimation of brain connectivity measured with EEG/MEG,” Computers in Biology and Medicine, in press. [11] T. Cover and J. Thomas, Elements of Information Theory. Wiley-Interscience, 2006. [12] M. Kaminski and K. Blinowska, “A new method of the description of the information flow in the brain structures,” Biological Cybernetics, vol. 65, no. 3, pp. 203–210, 1991. 154 [13] L. Baccala and K. Sameshima, “Partial directed coherence: a new concept in neural structure determination,” Biological Cybernetics, vol. 84, no. 6, pp. 463–474, 2001. [14] C. Porcaro, F. Zappasodi, P. Rossini, and F. Tecchio, “Choice of multivariate autoregressive model order affecting real network functional connectivity estimate,” Clinical Neurophysiology, vol. 120, no. 2, pp. 436–448, 2009. [15] E. Fanselow, K. Sameshima, L. Baccala, and M. Nicolelis, “Thalamic bursting in rats during different awake behavioral states,” Proceedings of the National Academy of Sciences, vol. 98, no. 26, pp. 15 330–15 335, 2001. [16] G. Supp, A. Schlögl, N. Trujillo-Barreto, M. Müller, and T. Gruber, “Directed cortical information flow during human object recognition: analyzing induced EEG gamma-band responses in brain’s source space,” PLoS ONE, vol. 2, no. 8, p. 684. [17] M. Kaminski, M. Ding, W. Truccolo, and S. Bressler, “Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of significance,” Biological Cybernetics, vol. 85, no. 2, pp. 145–157, 2001. [18] O. David, D. Cosmelli, and K. Friston, “Evaluation of different measures of functional connectivity using a neural mass model,” NeuroImage, vol. 21, no. 2, pp. 659–673, 2004. [19] R. Kus, M. Kaminski, and K. Blinowska, “Determination of EEG activity propagation: pair-wise versus multichannel estimate,” IEEE Transactions on Biomedical Engineering, vol. 51, no. 9, pp. 1501–1510, 2004. [20] A. Schlögl and G. Supp, “Analyzing event-related EEG data with multivariate autoregressive parameters,” Progress in brain research, vol. 159, no. 1, pp. 135–147, 2006. [21] M. Scherg, “Fundamentals of dipole source potential analysis,” Advances in Audiology, vol. 6, no. 1, pp. 40–69, 1990. [22] Z. Koles, “Trends in EEG source localization,” Electroencephalography and Clinical Neurophysiology, vol. 106, no. 2, pp. 127–137, 1998. [23] C. James and C. Hesse, “Independent component analysis for biomedical signals,” Physiological Measurement, vol. 26, no. 1, pp. 15–39, 2005. [24] T. Mima and M. Hallett, “Corticomuscular coherence: a review.” Journal of Clinical Neurophysiology, vol. 16, no. 6, pp. 501–511, 1999. [25] P. Grosse, M. Cassidy, and P. Brown, “EEG–EMG, MEG–EMG and EMG–EMG frequency analysis: physiological principles and clinical applications,” Clinical Neurophysiology, vol. 113, no. 10, pp. 1523–1531, 2002. [26] B. Conway, D. Halliday, S. Farmer, U. Shahani, P. Maas, A. Weir, and J. Rosenberg, “Synchronization between motor cortex and spinal motoneuronal pool during the performance of a maintained motor task in man,” The Journal of Physiology, vol. 489, no. 3, pp. 917–924, 1995. 155 [27] S. Salenius, K. Portin, M. Kajola, R. Salmelin, and R. Hari, “Cortical control of human motoneuron firing during isometric contraction,” Journal of Neurophysiology, vol. 77, no. 6, pp. 3401–3405, 1997. [28] D. Halliday, B. Conway, S. Farmer, and J. Rosenberg, “Using electroencephalography to study functional coupling between cortical activity and electromyograms during voluntary contractions in humans,” Neuroscience Letters, vol. 241, no. 1, pp. 5–8, 1998. [29] B. Feige, A. Aertsen, and R. Kristeva-Feige, “Dynamic synchronization between multiple cortical motor areas and muscle activity in phasic voluntary movements,” Journal of Neurophysiology, vol. 84, no. 5, pp. 2622–2629, 2000. [30] S. Salenius, S. Avikainen, S. Kaakkola, R. Hari, and P. Brown, “Defective cortical drive to muscle in Parkinson’s disease and its improvement with levodopa,” Brain, vol. 125, no. 3, pp. 491–500, 2002. [31] Y. Fang, J. Daly, J. Sun, K. Hvorat, E. Fredrickson, S. Pundik, V. Sahgal, and G. Yue, “Functional corticomuscular connection during reaching is weakened following stroke,” Clinical Neurophysiology, vol. 120, no. 5, pp. 994–1002, 2009. [32] V. Murthy and E. Fetz, “Coherent 25-to 35-Hz oscillations in the sensorimotor cortex of awake behaving monkeys,” Proceedings of the National Academy of Sciences of the United States of America, vol. 89, no. 12, pp. 5670–5674, 1992. [33] R. Bortel and P. Sovka, “EEG-EMG coherence enhancement,” Signal Processing, vol. 86, no. 7, pp. 1737–1751, 2006. [34] A. Samii, J. Nutt, and B. Ransom, “Parkinson’s disease,” The Lancet, vol. 363, no. 9423, pp. 1783–1793, 2004. [35] P. S. Canada, “Parkinson’s disease: Social and economic impact,” 2003. [36] Movement Disorder Society Task Force on Rating Scales for Parkinsons Disease, “The Unified Parkinson’s Disease Rating Scale (UPDRS): Status and recommendations,” Movement Disorders, vol. 18, no. 7, pp. 738–750, 2003. [37] A. Benabid, “Deep brain stimulation for parkinsons disease,” Current Opinion in Neurobiology, vol. 13, no. 6, pp. 696–706, 2003. [38] P. Brown, “Bad oscillations in Parkinson’s disease,” Journal of neural transmission. Supplementum, vol. 70, no. 1, pp. 27–30, 2006. [39] J. Dushanova, D. Philipova, and G. Nikolova, “Beta and gamma frequency-range abnormalities in parkinsonian patients under cognitive sensorimotor task,” Journal of the Neurological Sciences, vol. 293, no. 1-2, pp. 51–58, 2010. [40] A. Kühn, A. Tsui, T. Aziz, N. Ray, C. Brücke, A. Kupsch, G. Schneider, and P. Brown, “Pathological synchronisation in the subthalamic nucleus of patients with Parkinson’s disease relates to both bradykinesia and rigidity,” Experimental Neurology, vol. 215, no. 2, pp. 380–387, 2009. 156 [41] A. Kühn, F. Kempf, C. Brücke, L. Gaynor Doyle, I. Martinez-Torres, A. Pogosyan, T. Trottenberg, A. Kupsch, G. Schneider, M. Hariz, W. Vandenberghe, B. Nuttin, and P. Brown, “High-frequency stimulation of the subthalamic nucleus suppresses oscillatory beta activity in patients with Parkinson’s disease in parallel with improvement in motor performance,” The Journal of Neuroscience, vol. 28, no. 24, pp. 6165–6173, 2008. [42] P. Brown, “Bradykinesia and impairment of EEG desynchronization in Parkinson’s disease,” Movement Disorders, vol. 14, no. 3, pp. 423–429, 1999. [43] M. Asgari and I. Shafran, “Predicting severity of Parkinson’s disease from speech,” in Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC’10), Aug. 2010, pp. 5201–5204. [44] A. Tsanas, M. Little, P. McSharry, and L. Ramig, “Accurate telemonitoring of Parkinson’s disease progression by noninvasive speech tests,” IEEE Transactions on Biomedical Engineering, vol. 57, no. 4, pp. 884–893, 2010. [45] M. Ding, Y. Chen, and S. Bressler, Granger causality: Basic theory and application to neuroscience. Wiley, 2006, pp. 437–460. [46] P. Valdés-Sosa, J. Sánchez-Bornot, A. Lage-Castellanos, M. Vega-Hernández, J. Bosch-Bayard, L. Melie-Garcı́a, and E. Canales-Rodrı́guez, “Estimating brain functional connectivity with sparse multivariate autoregression,” Philosophical Transactions of the Royal Society B: Biological Sciences, vol. 360, no. 1457, pp. 969–981, 2005. [47] J. Fan and R. Li, “Variable selection via nonconcave penalized likelihood and its oracle properties,” Journal of the American Statistical Association, vol. 96, no. 456, pp. 1348–1360, 2001. [48] T. Jung, C. Humphries, T. Lee, S. Makeig, M. McKeown, V. Iragui, and T. Sejnowski, “Extended ICA removes artifacts from electroencephalographic recordings,” Advances in Neural Information Processing Systems, vol. 10, no. 1, pp. 894–900, 1998. [49] C. Başar-Eroglu, D. Strüber, P. Kruse, E. Başar, and M. Stadler, “Frontal gamma-band enhancement during multistable visual perception,” International Journal of Psychophysiology, vol. 24, no. 1-2, pp. 113–125, 1996. [50] P. Brown, “Oscillatory nature of human basal ganglia activity: relationship to the pathophysiology of Parkinson’s disease,” Movement Disorders, vol. 18, no. 4, pp. 357–363, 2003. [51] P. Brown, A. Oliviero, P. Mazzone, A. Insola, P. Tonali, and V. Di Lazzaro, “Dopamine dependency of oscillations between subthalamic nucleus and pallidum in Parkinson’s disease,” The Journal of Neuroscience, vol. 21, no. 3, pp. 1033–1038, 2001. 157 [52] R. Amirnovin, Z. Williams, G. Cosgrove, and E. Eskandar, “Visually guided movements suppress subthalamic oscillations in Parkinson’s disease patients,” The Journal of Neuroscience, vol. 24, no. 50, pp. 11 302–11 306, 2004. [53] A. Kühn, D. Williams, A. Kupsch, P. Limousin, M. Hariz, G. Schneider, K. Yarrow, and P. Brown, “Event-related beta desynchronization in human subthalamic nucleus correlates with motor performance,” Brain, vol. 127, no. 4, pp. 735–746, 2004. [54] D. Williams, A. Kuhn, A. Kupsch, M. Tijssen, G. van Bruggen, H. Speelman, G. Hotton, K. Yarrow, and P. Brown, “Behavioural cues are associated with modulations of synchronous oscillations in the human subthalamic nucleus,” Brain, vol. 126, no. 9, pp. 1975–1985, 2003. [55] R. Levy, P. Ashby, W. Hutchison, A. Lang, A. Lozano, and J. Dostrovsky, “Dependence of subthalamic nucleus oscillations on movement and dopamine in Parkinson’s disease,” Brain, vol. 125, no. 6, pp. 1196–1209, 2002. [56] D. Williams, M. Tijssen, G. van Bruggen, A. Bosch, A. Insola, V. DiLazzaro, P. Mazzone, A. Oliviero, A. Quartarone, H. Speelman, and P. Brown, “Dopamine-dependent changes in the functional connectivity between basal ganglia and cerebral cortex in humans,” Brain, vol. 125, no. 7, pp. 1558–1569, 2002. [57] C. Chen, V. Litvak, T. Gilbertson, A. Kühn, C. Lu, S. Lee, C. Tsai, S. Tisch, P. Limousin, M. Hariz, and P. Brown, “Excessive synchronization of basal ganglia neurons at 20 Hz slows movement in Parkinson’s disease,” Experimental Neurology, vol. 205, no. 1, pp. 214–221, 2007. [58] G. Arbuthnott and M. Garcia-Munoz, “Dealing with the devil in the detail - some thoughts about the next model of the basal ganglia,” Parkinsonism & Related Disorders, vol. Suppl 3, no. 1, pp. 139–142, 2009. [59] P. Brown and D. Williams, “Basal ganglia local field potential activity: character and functional significance in the human,” Clinical Neurophysiology, vol. 116, no. 11, pp. 2510–2519, 2005. [60] M. Cassidy, P. Mazzone, A. Oliviero, A. Insola, P. Tonali, V. Di Lazzaro, and P. Brown, “Movement-related changes in synchronization in the human basal ganglia,” Brain, vol. 125, no. 6, pp. 1235–1246, 2002. [61] P. Gatev, O. Darbin, and T. Wichmann, “Oscillations in the basal ganglia under normal conditions and in movement disorders,” Movement Disorders, vol. 21, no. 10, pp. 1566–1577, 2006. [62] S. Li, G. Arbuthnott, M. Jutras, J. Goldberg, and D. Jaeger, “Resonant antidromic cortical circuit activation as a consequence of high-frequency subthalamic deep-brain stimulation,” Journal of Neurophysiology, vol. 98, no. 6, pp. 3525–3537, 2007. [63] C. Dejean, B. Hyland, and G. Arbuthnott, “Cortical effects of subthalamic stimulation correlate with behavioral recovery from dopamine antagonist induced akinesia,” Cerebral Cortex, vol. 19, no. 5, pp. 1055–1063, 2009. 158 [64] V. Gradinaru, M. Mogri, K. Thompson, J. Henderson, and K. Deisseroth, “Optical deconstruction of parkinsonian neural circuitry,” Science, vol. 324, no. 5925, pp. 354–359, 2009. [65] L. Doyle Gaynor, A. Kühn, M. Dileone, V. Litvak, A. Eusebio, A. Pogosyan, A. Androulidakis, S. Tisch, P. Limousin, A. Insola, P. Mazzone, V. Di Lazzaro, and P. Brown, “Suppression of beta oscillations in the subthalamic nucleus following cortical stimulation in humans,” European Journal of Neuroscience, vol. 28, no. 8, pp. 1686–1695, 2008. [66] N. Fogelson, D. Williams, M. Tijssen, G. van Bruggen, H. Speelman, and P. Brown, “Different functional loops between cerebral cortex and the subthalmic area in Parkinson’s disease,” Cerebral Cortex, vol. 16, no. 1, pp. 64–75, 2006. [67] F. Klostermann, V. Nikulin, A. Kühn, F. Marzinzik, M. Wahl, A. Pogosyan, A. Kupsch, G. Schneider, P. Brown, and G. Curio, “Task-related differential dynamics of EEG alpha-and beta-band synchronization in cortico-basal motor structures,” European Journal of Neuroscience, vol. 25, no. 5, pp. 1604–1615, 2007. [68] E. Lalo, S. Thobois, A. Sharott, G. Polo, P. Mertens, A. Pogosyan, and P. Brown, “Patterns of bidirectional communication between cortex and basal ganglia during movement in patients with Parkinson disease,” The Journal of Neuroscience, vol. 28, no. 12, pp. 3008–3016, 2008. [69] J. Marsden, P. Limousin-Dowsey, P. Ashby, P. Pollak, and P. Brown, “Subthalamic nucleus, sensorimotor cortex and muscle interrelationships in Parkinson’s disease,” Brain, vol. 124, no. 2, pp. 378–388, 2001. [70] M. Cassidy and P. Brown, “Task-related EEG-EEG coherence depends on dopaminergic activity in Parkinson’s disease,” Neuroreport, vol. 12, no. 4, pp. 703–707, 2001. [71] P. Silberstein, A. Pogosyan, A. Kühn, G. Hotton, S. Tisch, A. Kupsch, P. Dowsey-Limousin, M. Hariz, and P. Brown, “Cortico-cortical coupling in Parkinson’s disease and its modulation by therapy,” Brain, vol. 128, no. 6, pp. 1277–1291, 2005. [72] L. Astolfi, F. Cincotti, D. Mattia, M. Marciani, L. Baccala, F. de Vico Fallani, S. Salinari, M. Ursino, M. Zavaglia, L. Ding, J. Edgar, G. Miller, B. He, and F. Babiloni, “Comparison of different cortical connectivity estimators for high-resolution EEG recordings,” Human Brain Mapping, vol. 28, no. 2, pp. 143–157, 2007. [73] H. Witte, M. Ungureanu, C. Ligges, D. Hemmelmann, T. Wstenberg, J. Reichenbach, L. Astolfi, F. Babiloni, and L. Leistritz, “Signal informatics as an advanced integrative concept in the framework of medical informatics. New trends demonstrated by examples derived from neuroscience,” Methods of Information in Medicine, vol. 48, no. 1, pp. 18–28, 2009. 159 [74] Y. Sun, H. Zhang, T. Feng, Y. Qiu, Y. Zhu, and S. Tong, “Early cortical connective network relating to audiovisual stimulation by partial directed coherence analysis,” IEEE Transactions on Biomedical Engineering, vol. 56, no. 11, pp. 2721–2724, 2009. [75] P. Nunez, R. Srinivasan, A. Westdorp, R. Wijesinghe, D. Tucker, R. Silberstein, and P. Cadusch, “EEG coherency I: statistics, reference electrode, volume conduction, Laplacians, cortical imaging, and interpretation at multiple scales,” Electroencephalography and Clinical Neurophysiology, vol. 103, no. 5, pp. 499–515, 1997. [76] R. Srinivasan, W. Winter, J. Ding, and P. Nunez, “EEG and MEG coherence: measures of functional connectivity at distinct spatial scales of neocortical dynamics,” Journal of Neuroscience Methods, vol. 166, no. 1, pp. 41–52, 2007. [77] W. Winter, P. Nunez, J. Ding, and R. Srinivasan, “Comparison of the effect of volume conduction on EEG coherence with the effect of field spread on MEG coherence,” Statistics in Medicine, vol. 26, no. 1, pp. 3946–3957, 2007. [78] S. Palmer, L. Eigenraam, T. Hoque, R. McCaig, A. Troiano, and M. McKeown, “Levodopa-sensitive, dynamic changes in effective connectivity during simultaneous movements in Parkinson’s disease,” Neuroscience, vol. 158, no. 2, pp. 693–704, 2009. [79] S. Palmer, B. Ng, R. Abugharbieh, L. Eigenraam, and M. McKeown, “Motor reserve and novel area recruitment: amplitude and spatial characteristics of compensation in Parkinson’s disease,” European Journal of Neuroscience, vol. 29, no. 11, pp. 2187–2196, 2009. [80] S. Palmer, J. Li, Z. Wang, and M. McKeown, “Joint amplitude and connectivity compensatory mechanisms in Parkinson’s disease,” Neuroscience, vol. 166, no. 4, pp. 1110–1118, 2010. [81] D. Brainard, “The psychophysics toolbox,” Spatial Vision, vol. 10, no. 4, pp. 433–436, 1997. [82] G. Gómez-Herrero, “Automatic Artifact Removal (AAR) toolbox v1. 3 (Release 09.12. 2007) for MATLAB,” 2007, http://www.cs.tut.fi/∼gomezher/projects/eeg/aar.htm. [83] 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. [84] R. Homan, J. Herman, and P. Purdy, “Cerebral location of international 10-20 system electrode placement,” Electroencephalography and Clinical Neurophysiology, vol. 66, no. 4, pp. 376–382, 1987. [85] H. Steinmetz, G. Furst, and B. Meyer, “Craniocerebral topography within the international 10-20 system,” Electroencephalography and Clinical Neurophysiology, vol. 72, no. 6, pp. 499–506, 1989. 160 [86] L. Doyle, K. Yarrow, and P. Brown, “Lateralization of event-related beta desynchronization in the EEG during pre-cued reaction time tasks,” Clinical Neurophysiology, vol. 116, no. 8, pp. 1879–1888, 2005. [87] C. Kranczioch, S. Athanassiou, S. Shen, G. Gao, and A. Sterr, “Short-term learning of a visually guided power-grip task is associated with dynamic changes in EEG oscillatory activity,” Clinical Neurophysiology, vol. 119, no. 6, pp. 1419–1430, 2008. [88] P. Praamstra and P. Pope, “Slow brain potential and oscillatory EEG manifestations of impaired temporal preparation in Parkinson’s disease,” Journal of Neurophysiology, vol. 98, no. 5, pp. 2848–2857, 2007. [89] D. Stoffers, J. Bosboom, E. Wolters, C. Stam, and H. Berendse, “Dopaminergic modulation of cortico-cortical functional connectivity in Parkinson’s disease: An MEG study,” Experimental Neurology, vol. 213, no. 1, pp. 191–195, 2008. [90] M. Moazami-Goudarzi, J. Sarnthein, L. Michels, R. Moukhtieva, and D. Jeanmonod, “Enhanced frontal low and high frequency power and synchronization in the resting EEG of parkinsonian patients,” NeuroImage, vol. 41, no. 3, pp. 985–997, 2008. [91] N. Crone, A. Sinai, and A. Korzeniewska, “High-frequency gamma oscillations and human brain mapping with electrocorticography,” Progress in Brain Research, vol. 159, no. 1, pp. 275–295, 2006. [92] I. Bodis-Wollner, “Visual deficits related to dopamine deficiency in experimental animals and Parkinson’s disease patients,” Trends in Neurosciences, vol. 13, no. 7, pp. 296–302, 1990. [93] R. Rodnitzky, “Visual dysfunction in Parkinson’s disease,” Clinical Neuroscience, vol. 5, no. 2, pp. 102–106, 1998. [94] M. Silva, P. Faria, F. Regateiro, V. Forjaz, C. Januario, A. Freire, and M. Castelo-Branco, “Independent patterns of damage within magno-, parvo- and koniocellular pathways in Parkinson’s disease,” Brain, vol. 128, no. 10, pp. 2260–2271, 2005. [95] A. Antal, F. Bandini, S. Keri, and I. Bodis-Wollner, “Visuo-cognitive dysfunctions in Parkinson’s disease,” Clinical Neuroscience, vol. 5, no. 2, pp. 147–152, 1998. [96] Y. Abe, T. Kachi, T. Kato, Y. Arahata, T. Yamada, Y. Washimi, K. Iwai, K. Ito, N. Yanagisawa, and G. Sobue, “Occipital hypoperfusion in Parkinson’s disease without dementia: correlation to impaired cortical visual processing,” Journal of Neurology, Neurosurgery & Psychiatry, vol. 74, no. 4, pp. 419–422, 2003. [97] J. Pereira, C. Junque, M. Marti, B. Ramirez-Ruiz, N. Bargallo, and E. Tolosa, “Neuroanatomical substrate of visuospatial and visuoperceptual impairment in Parkinson’s disease,” Movement Disorders, vol. 24, no. 8, pp. 1193–1199, 2009. 161 [98] M. Jahanshahi, I. Jenkins, R. Brown, C. Marsden, R. Passingham, and D. Brooks, “Self-initiated versus externally triggered movements: I. An investigation using measurement of regional cerebral blood flow with PET and movement-related potentials in normal and Parkinson’s disease subjects,” Brain, vol. 118, no. 4, pp. 913–933, 1995. [99] M. Lewis, C. Slagle, A. Smith, Y. Truong, P. Bai, M. McKeown, R. Mailman, A. Belger, and X. Huang, “Task specific influences of Parkinson’s disease on the striato-thalamo-cortical and cerebello-thalamo-cortical motor circuitries,” Neuroscience, vol. 147, no. 1, pp. 224–235, 2007. [100] R. Siegert, D. Harper, F. Cameron, and D. Abernethy, “Self-initiated versus externally cued reaction times in Parkinson’s disease,” Journal of Clinical and Experimental Neuropsychology, vol. 24, no. 2, pp. 146–153, 2002. [101] K. Werheid, I. Koch, K. Reichert, and M. Brass, “Impaired self-initiated task preparation during task switching in Parkinson’s disease,” Neuropsychologia, vol. 45, no. 2, pp. 273–281, 2007. [102] S. Sen, A. Kawaguchi, Y. Truong, M. Lewis, and X. Huang, “Dynamic changes in cerebello-thalamo-cortical motor circuitry during progression of Parkinson’s disease,” Neuroscience, vol. 166, no. 2, pp. 712–719, 2010. [103] T. Taniwaki, A. Okayama, T. Yoshiura, O. Togao, Y. Nakamura, T. Yamasaki, K. Ogata, H. Shigeto, Y. Ohyagi, J. Kira, and S. Tobimatsu, “Functional network of the basal ganglia and cerebellar motor loops in vivo: Different activation patterns between self-initiated and externally triggered movements,” NeuroImage, vol. 31, no. 2, pp. 745–753, 2006. [104] P. Damier, “Drug-induced dyskinesias,” Current Opinion in Neurology, vol. 22, no. 4, pp. 394–399, 2009. [105] A. Schapira, M. Emre, P. Jenner, and W. Poewe, “Levodopa in the treatment of Parkinson’s disease,” European Journal of Neurology, vol. 16, no. 9, pp. 982–989, 2009. [106] R. Ferri, F. Rundo, O. Bruni, M. Terzano, and C. Stam, “The functional connectivity of different EEG bands moves towards small-world network organization during sleep,” Clinical Neurophysiology, vol. 119, no. 9, pp. 2026–2036, 2008. [107] S. Micheloyannis, E. Pachou, C. Stam, M. Breakspear, P. Bitsios, M. Vourkas, S. Erimaki, and M. Zervakis, “Small-world networks and disturbed functional connectivity in schizophrenia,” Schizophrenia Research, vol. 87, no. 1-3, pp. 60–66, 2006. [108] J. Li, Z. Wang, S. Palmer, and M. McKeown, “Dynamic Bayesian network modeling of fMRI: a comparison of group-analysis methods,” NeuroImage, vol. 41, no. 2, pp. 398–407, 2008. 162 [109] G. Tropini, J. Chiang, Z. J. Wang, E. Ty, and M. J. McKeown, “Altered directional connectivity in Parkinson’s disease during performance of a visually guided task,” NeuroImage, vol. 56, no. 4, pp. 2144–2156, Jun. 2011. [110] J. Bosboom, D. Stoffers, C. Stam, B. van Dijk, J. Verbunt, H. Berendse, and E. Wolters, “Resting state oscillatory brain dynamics in Parkinson’s disease: An MEG study,” Clinical Neurophysiology, vol. 117, no. 11, pp. 2521–2531, 2006. [111] D. Stoffers, J. Bosboom, J. Deijen, E. Wolters, H. Berendse, and C. Stam, “Slowing of oscillatory brain activity is a stable characteristic of Parkinson’s disease without dementia,” Brain, vol. 130, no. 7, pp. 1847–1860, 2007. [112] R. Tibshirani, “Regression shrinkage and selection via the lasso,” Journal of the Royal Statistical Society. Series B (Methodological), vol. 58, no. 1, pp. 267–288, 1996. [113] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, “Least angle regression,” The Annals of statistics, vol. 32, no. 2, pp. 407–499, 2004. [114] A. D.E. and C. K.E., “Vestibular system: the many facets of a multimodal sense,” Annual Review of Neuroscience, vol. 31, no. 1, pp. 125–150, 2008. [115] W. Guldin and O. Grüsser, “Is there a vestibular cortex?” Trends in Neurosciences, vol. 21, no. 6, pp. 254–259, 1998. [116] C. Lopez and O. Blanke, “The thalamocortical vestibular system in animals and humans,” Brain Research Reviews, vol. 67, no. 1, pp. 119–146, 2011. [117] R. Fitzpatrick and B. Day, “Probing the human vestibular system with galvanic stimulation,” Journal of Applied Physiology, vol. 96, no. 6, pp. 2301–2316, 2004. [118] D. Wilkinson, S. Nicholls, C. Pattenden, P. Kilduff, and W. Milberg, “Galvanic vestibular stimulation speeds visual memory recall,” Experimental Brain Research, vol. 189, no. 2, pp. 243–248, 2008. [119] B. Day, “Galvanic vestibular stimulation: new uses for an old tool,” The Journal of Physiology, vol. 517, no. 3, p. 631, 1999. [120] Y. Yamamoto, Z. Struzik, R. Soma, K. Ohashi, and S. Kwak, “Noisy vestibular stimulation improves autonomic and motor responsiveness in central neurodegenerative disorders,” Annals of Neurology, vol. 58, no. 2, pp. 175–181, 2005. [121] W. Pan, R. Soma, S. Kwak, and Y. Yamamoto, “Improvement of motor functions by noisy vestibular stimulation in central neurodegenerative disorders,” Journal of Neurology, vol. 255, no. 11, pp. 1657–1661, 2008. [122] C. Hamani, J. Saint-Cyr, J. Fraser, M. Kaplitt, and A. Lozano, “The subthalamic nucleus in the context of movement disorders,” Brain, vol. 127, no. 1, pp. 4–20, 2004. 163 [123] R. Llinás, U. Ribary, D. Jeanmonod, E. Kronberg, and P. Mitra, “Thalamocortical dysrhythmia: a neurological and neuropsychiatric syndrome characterized by magnetoencephalography,” Proceedings of the National Academy of Sciences, vol. 96, no. 26, pp. 15 222–15 227, 1999. [124] T. Jung, S. Makeig, C. Humphries, T. Lee, M. Mckeown, V. Iragui, and T. Sejnowski, “Removing electroencephalographic artifacts by blind source separation,” Psychophysiology, vol. 37, no. 2, pp. 163–178, 2000. [125] P. Rousseeuw and A. Leroy, Robust Regression and Outlier Detection, 1st ed. New York, NY: John Wiley & Sons, Inc., 1987. [126] M. Yuan and Y. Lin, “Model selection and estimation in regression with grouped variables,” Journal of the Royal Statistical Society: Series B, vol. 68, no. 1, pp. 49–67, 2006. [127] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 1.21,” 2010. [128] M. Nuwer, G. Comi, R. Emerson, A. Fuglsang-Frederiksen, J. Guérit, H. Hinrichs, A. Ikeda, F. Luccas, and P. Rappelsburger, “IFCN standards for digital recording of clinical EEG,” Electroencephalography and Clinical Neurophysiology, vol. 106, no. 3, pp. 259–261, 1998. [129] S. Balter, M. Castelijns, R. Stokroos, and H. Kingma, “Galvanic-induced body sway in vestibular schwannoma patients: evidence for stimulation of the central vestibular system,” Acta Otolaryngol, vol. 124, no. 9, pp. 1015–1021, 2004. [130] H. Meng, P. May, J. Dickman, and D. Angelaki, “Vestibular signals in primate thalamus: properties and origins,” The Journal of Neuroscience, vol. 27, no. 50, pp. 13 590–13 602, 2007. [131] G. Percheron, C. Francois, B. Talbi, J. Yelnik, and G. Fenelon, “The primate motor thalamus,” Brain Research Reviews, vol. 22, no. 2, pp. 93–181, 1996. [132] S. Sherman, “Thalamic relays and cortical functioning,” Progress in Brain Research, vol. 149, no. 1, pp. 107–126, 2005. [133] M. Emri, M. Kisely, Z. Lengyel, L. Balkay, T. Márián, L. Mikó, E. Berényi, I. Sziklai, L. Trón, and Á. Tóth, “Cortical projection of peripheral vestibular signaling,” Journal of Neurophysiology, vol. 89, no. 5, pp. 2639–2646, 2003. [134] O. Fasold, M. von Brevern, M. Kuhberg, C. Ploner, A. Villringer, T. Lempert, and R. Wenzel, “Human vestibular cortex as identified with caloric stimulation in functional magnetic resonance imaging,” NeuroImage, vol. 17, no. 3, pp. 1384–1393, 2002. [135] E. Lobel, J. Kleine, A. Leroy-Willig, P. Moortele, D. Bihan, O. Grüsser, and A. Berthoz, “Cortical areas activated by bilateral galvanic vestibular stimulation,” Annals of the New York Academy of Sciences, vol. 871, no. 1, pp. 313–323, 1999. 164 [136] E. Vitte, C. Derosier, Y. Caritu, A. Berthoz, D. Hasboun, and D. Soulié, “Activation of the hippocampal formation by vestibular stimulation: a functional magnetic resonance imaging study,” Experimental Brain Research, vol. 112, no. 3, pp. 523–526, 1996. [137] M. Kluge, S. Beyenburg, G. Fernández, and C. Elger, “Epileptic vertigo: evidence for vestibular representation in human frontal cortex,” Neurology, vol. 55, no. 12, pp. 1906–1908, 2000. [138] A. Cauquil and B. Day, “Galvanic vestibular stimulation modulates voluntary movement of the human upper body,” The Journal of Physiology, vol. 513, no. 2, pp. 611–619, 1998. [139] B. Day, A. Severac Cauquil, L. Bartolomei, M. Pastor, and I. Lyon, “Human body-segment tilts induced by galvanic stimulation: a vestibularly driven balance protection mechanism,” The Journal of Physiology, vol. 500, no. 3, pp. 661–672, 1997. [140] L. Bent, B. McFadyen, V. Merkley, P. Kennedy, and J. Inglis, “Magnitude effects of galvanic vestibular stimulation on the trajectory of human gait,” Neuroscience Letters, vol. 279, no. 3, pp. 157–160, 2000. [141] R. Fitzpatrick, D. Wardman, and J. Taylor, “Effects of galvanic vestibular stimulation during human walking,” The Journal of Physiology, vol. 517, no. 3, pp. 931–939, 1999. [142] A. Eusebio and P. Brown, “Synchronisation in the beta frequency-band–The bad boy of parkinsonism or an innocent bystander?” Experimental Neurology, vol. 217, no. 1, pp. 1–3, 2009. [143] M. McDonnell and L. Ward, “The benefits of noise in neural systems: bridging theory and experiment,” Nature Reviews Neuroscience, vol. 12, no. 7, pp. 415–426, 2011. [144] F. Moss, L. Ward, and W. Sannita, “Stochastic resonance and sensory information processing: a tutorial and review of application,” Clinical Neurophysiology, vol. 115, no. 2, pp. 267–281, 2004. [145] B. Cohen, G. Martinelli, D. Ogorodnikov, Y. Xiang, T. Raphan, G. Holstein, and S. Yakushin, “Sinusoidal galvanic vestibular stimulation (sGVS) induces a vasovagal response in the rat,” Experimental Brain Research, vol. 210, no. 1, pp. 45–55, 2011. [146] P. Smith, L. Geddes, J. Baek, C. Darlington, and Y. Zheng, “Modulation of memory by vestibular lesions and galvanic vestibular stimulation,” Frontiers in Neurology, vol. 1, pp. 1–8, 2010. [147] F. Babiloni, C. Babiloni, F. Carducci, L. Fattorini, P. Onorati, and A. Urbano, “Spline Laplacian estimate of EEG potentials over a realistic magnetic resonance-constructed scalp surface model,” Electroencephalography and Clinical Neurophysiology, vol. 98, no. 4, pp. 363–373, 1996. 165 [148] G. Gómez-Herrero, M. Atienza, K. Egiazarian, and J. Cantero, “Measuring directional coupling between EEG sources,” NeuroImage, vol. 43, no. 3, pp. 497–508, 2008. [149] S. Haufe, R. Tomioka, G. Nolte, K. Müller, and M. Kawanabe, “Modeling sparse connectivity between underlying brain sources for EEG/MEG,” IEEE Transactions on Biomedical Engineering, vol. 57, no. 8, pp. 1954–1963, 2010. [150] S. Baillet, J. Mosher, and R. Leahy, “Electromagnetic brain mapping,” IEEE Signal Processing Magazine, vol. 18, no. 6, pp. 14–30, 2001. [151] L. Zhukov, D. Weinstein, and C. Johnson, “Independent component analysis for EEG source localization,” IEEE Engineering in Medicine and Biology Magazine, vol. 19, no. 3, pp. 87–96, 2002. [152] A. Hyvärinen, P. Ramkumar, L. Parkkonen, and R. Hari, “Independent component analysis of short-time Fourier transforms for spontaneous EEG/MEG analysis,” NeuroImage, vol. 49, no. 1, pp. 257–271, 2010. [153] A. Hyvärinen, K. Zhang, S. Shimizu, and P. Hoyer, “Estimation of a structural vector autoregression model using non-Gaussianity,” The Journal of Machine Learning Research, vol. 11, no. 1, pp. 1709–1731, 2010. [154] B. Cheung, B. Riedner, G. Tononi, and B. Van Veen, “Estimation of cortical connectivity from EEG using state-space models,” IEEE Transactions on Biomedical Engineering, vol. 57, no. 9, pp. 2122–2134, 2010. [155] W. Penny, R. Everson, and S. Roberts, “Hidden Markov independent component analysis,” Advances in Independent Component Analysis, pp. 3–22, 2000. [156] S. Choi, A. Cichocki, and S. Amari, “Flexible independent component analysis,” The Journal of VLSI Signal Processing, vol. 26, no. 1, pp. 25–38, 2000. [157] G. Casella, R. Berger, and R. Berger, Statistical Inference, 2nd ed. Duxbury Press, 2002. [158] J. Cantero, M. Atienza, G. Gomez-Herrero, A. Cruz-Vadell, E. Gil-Neciga, R. Rodriguez-Romero, and D. Garcia-Solis, “Functional integrity of thalamocortical circuits differentiates normal aging from mild cognitive impairment,” Human Brain Mapping, vol. 30, no. 12, pp. 3944–3957, 2009. [159] J. Himberg, A. Hyvärinen, and F. Esposito, “Validating the independent components of neuroimaging time series via clustering and visualization,” NeuroImage, vol. 22, no. 3, pp. 1214–1222, 2004. [160] G. Gómez-Herrero, “MVAR-ICA toolbox for MATLAB,” 2008, http://www.cs.tut.fi/∼gomezher/projects/eeg/mvarica.htm. [161] D. Brody, F. Terry, and R. Ideker, “Eccentric dipole in a spherical medium: generalized expression for surface potentials,” IEEE Transactions on Biomedical Engineering, vol. 20, no. 2, pp. 141–143, 1973. 166 [162] P. Geladi and B. Kowalski, “Partial least-squares regression: a tutorial,” Analytica Chimica Acta, vol. 185, no. 1, pp. 1–17, 1986. [163] P. Tichavsky and Z. Koldovsky, “Optimal pairing of signal components separated by blind techniques,” IEEE Signal Processing Letters, vol. 11, no. 2, pp. 119–122, 2004. [164] T. Schneider and A. Neumaier, “Algorithm 808: ARfitA Matlab package for the estimation of parameters and eigenmodes of multivariate autoregressive models,” ACM Transactions on Mathematical Software (TOMS), vol. 27, no. 1, pp. 58–65, 2001. [165] P. Wang, C. Hung, S. Kwan, S. Teng, and B. Soong, “Early detection of periodic sharp wave complexes on EEG by independent component analysis in patients with Creutzfeldt-Jakob disease,” Journal of Clinical Neurophysiology, vol. 25, no. 1, pp. 25–31, 2008. [166] M. Rivlin-Etzion, O. Marmor, G. Heimer, A. Raz, A. Nini, and H. Bergman, “Basal ganglia oscillations and pathophysiology of movement disorders,” Current Opinion in Neurobiology, vol. 16, no. 6, pp. 629–637, 2006. [167] T. Jung, S. Makeig, M. McKeown, A. Bell, T. Lee, and T. Sejnowski, “Imaging brain dynamics using independent component analysis,” Proceedings of the IEEE, vol. 89, no. 7, pp. 1107–1122, 2001. [168] R. Vigário, J. Sarela, V. Jousmiki, M. Hamalainen, and E. Oja, “Independent component approach to the analysis of EEG and MEG recordings,” IEEE Transactions on Biomedical Engineering, vol. 47, no. 5, pp. 589–593, 2000. [169] S. Makeig, M. Westerfield, T. Jung, S. Enghoff, J. Townsend, E. Courchesne, and T. Sejnowski, “Dynamic brain sources of visual evoked responses,” Science, vol. 295, no. 5555, pp. 690–694, 2002. [170] J. Ylipaavalniemi and R. Vigário, “Analyzing consistency of independent components: An fMRI illustration,” NeuroImage, vol. 39, no. 1, pp. 169–180, 2008. [171] E. Learned-Miller and J. Fisher, III, “ICA using spacings estimates of entropy,” The Journal of Machine Learning Research, vol. 4, no. 7-8, pp. 1271–1295, 2003. [172] H. Stögbauer, A. Kraskov, S. Astakhov, and P. Grassberger, “Least-dependent-component analysis based on mutual information,” Physical Review E, vol. 70, no. 6, pp. 1–17, 2004. [173] H. Shen, S. Jegelka, and A. Gretton, “Fast kernel-based independent component analysis,” IEEE Transactions on Signal Processing, vol. 57, no. 9, pp. 3498–3511, 2009. [174] S. Baker, E. Olivier, and R. Lemon, “Coherent oscillations in monkey motor cortex and hand muscle EMG show task-dependent modulation,” The Journal of Physiology, vol. 501, no. 1, pp. 225–241, 1997. 167 [175] A. McIntosh, F. Bookstein, J. Haxby, and C. Grady, “Spatial pattern analysis of functional brain images using partial least squares,” NeuroImage, vol. 3, no. 3, pp. 143–157, 1996. [176] E. Martinez-Montes, P. Valdés-Sosa, F. Miwakeichi, R. Goldman, and M. Cohen, “Concurrent EEG/fMRI analysis by multiway partial least squares,” NeuroImage, vol. 22, no. 3, pp. 1023–1034, 2004. [177] R. Bro, “Multiway calibration. Multilinear PLS.” Journal of Chemometrics, vol. 10, pp. 47–61, 1996. [178] S. Wold, N. Kettaneh, and K. Tjessem, “Hierarchical multiblock PLS and PC models for easier model interpretation and as an alternative to variable selection,” Journal of Chemometrics, vol. 10, no. 5-6, pp. 463–482, 1998. [179] E. Clancy, E. Morin, and R. Merletti, “Sampling, noise-reduction and amplitude estimation issues in surface electromyography,” Journal of Electromyography and Kinesiology, vol. 12, no. 1, pp. 1–16, 2002. [180] P. Good, Permutation Tests: A Practical Guide to Resampling Methods for Testing Hypotheses, 2nd ed., ser. Springer Series in Statistics. Springer-Verlag, Heidelberg, 2000. [181] J. Rosenberg, A. Amjad, P. Breeze, D. Brillinger, and D. Halliday, “The Fourier approach to the identification of functional coupling between neuronal spike trains.” Progress in Biophysics and Molecular Biology, vol. 53, no. 1, pp. 1–31, 1989. [182] P. Praamstra, D. Stegeman, A. Cools, and M. Horstink, “Reliance on external cues for movement initiation in Parkinson’s disease. Evidence from movement-related potentials.” Brain, vol. 121, no. 1, pp. 167–177, 1998. [183] J. K. R. Stevenson, M. M. K. Oishi, S. Farajian, E. Cretu, E. Ty, and M. J. McKeown, “Response to sensory uncertainty in Parkinsons disease: a marker of cerebellar dysfunction?” European Journal of Neuroscience, vol. 33, no. 2, pp. 298–305, 2011. [184] R. Henson, D. Wakeman, V. Litvak, and K. Friston, “A parametric empirical Bayesian framework for the EEG/MEG inverse problem: generative models for multi-subject and multi-modal integration,” Frontiers in Human Neuroscience, vol. 5, no. 76, pp. 1–14, 2011. [185] V. Calhoun, T. Adali, M. Stevens, K. Kiehl, and J. Pekar, “Semi-blind ICA of fMRI : a method for utilizing hypothesis-derived time courses in a spatial ICA analysis,” NeuroImage, vol. 25, pp. 527–538, 2005. 168 Appendix A Estimation of Multiblock PLS Algorithm 2 Multiblock PLS algorithm 1: k← 1 2: Choose start uk, e.g., uk← the first principal component of some response blockY b. 3: repeat 4: for b= 1 to B do 5: wb,k←X Tbuk/ ‖ uTk uk ‖, normalize wb,k to ‖wb,k ‖= 1 6: rb,k←X bwb,k/wTb,kwb,k 7: end for 8: Rk← [ r1,k,r2,k, . . . ,rB,k ] 9: wsup←RTk uk/ ‖ uTk uk ‖, normalize wsup to ‖wsup ‖= 1 10: t k←Rkwsup/wTsupwsup 11: for b= 1 to B do 12: qb,k←Y Tb t k/ ‖ tTk t k ‖, normalize qb,k to ‖ qb,k ‖= 1 13: sb,k←Y bqb,k/qTb,kqb,k 14: end for 15: Sk← [ s1,k,s2,k, . . . ,sB,k ] 16: qsup← STk t k/ ‖ tTk t k ‖, normalize qsup to ‖ qsup ‖= 1 17: uk← SkqTsup/qTsupqsup 18: until t k converges 19: for b= 1 to B do 20: pb,k←X Tb t k/ ‖ tTk t k ‖ 21: X b←X b−rb,kpTb,k 22: Y b←Y b−sb,kqTb,k 23: end for 24: k← k+1. Repeat from Step 2 until X b andY b are properly described. 169
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Cortico-cortical and cortico-muscular connectivity...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Cortico-cortical and cortico-muscular connectivity analysis : methods and application to Parkinson's… Chiang, Joyce Hsien-yin 2012
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Cortico-cortical and cortico-muscular connectivity analysis : methods and application to Parkinson's disease |
Creator |
Chiang, Joyce Hsien-yin |
Publisher | University of British Columbia |
Date Issued | 2012 |
Description | The concept of brain connectivity provides a new perspective to the understanding of the mechanism underlying brain functions, complementing the traditional approach of analyzing neural activity of isolated regions. Among the existing connectivity analysis techniques, multivariate autoregressive (mAR)-based measures are of great interest for their ability to characterize both directionality and spectral property of cortical interactions. Yet, the direct estimation of mAR-based connectivity from scalp electroencephalogram (EEG) is confounded by volume conduction, statistical instability and inter-subject variability. In this thesis, we propose novel signal processing methods to enhance the existing mAR-based connectivity methods. First, we explore incorporating sparsity constraints into the mAR formulation at both subject level and group level using LASSO-based regression. We show by simulation that sparse mAR yields more stable and accurate connectivity estimates compared to the traditional, non-sparse approach. Furthermore, the group-wise sparsity simplifies the inference of group-level connectivity patterns from multi-subject data. To mitigate the effect of volume conduction, we investigate source-level connectivity and propose a state-space generalized mAR framework to jointly model the mixing effect of volume conduction and causal relationships between underlying neural sources. By jointly estimating the mixing process and mAR model parameters, the proposed technique demonstrates improved connectivity estimation performance. Finally, we expanded our connectivity analysis to cortico-muscular level by modeling the relationships between EEG and simultaneously recorded electromyography (EMG) data using a multiblock partial least square (mbPLS) framework. The hierarchical construction of the mbPLS framework provides a natural way to model multi-subject, multi-modal data, enabling the identification of maximally covarying common patterns from EEG and EMG across subjects. Applications of the proposed techniques to EEG and EMG data of healthy and Parkinson's disease (PD) subjects demonstrate that directional connectivity analysis is a more sensitive technique than traditional univariate spectral analysis in revealing complex effects of motor tasks and disease. Moreover, alternations in connectivity accurately predict disease severity in PD. These new analysis tools allow a better understanding of brain function and provide a basis for developing objective measures to assess progression of neurological diseases. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2012-04-19 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0103456 |
URI | http://hdl.handle.net/2429/42126 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2012-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2012_spring_chiang_joyce.pdf [ 6.92MB ]
- Metadata
- JSON: 24-1.0103456.json
- JSON-LD: 24-1.0103456-ld.json
- RDF/XML (Pretty): 24-1.0103456-rdf.xml
- RDF/JSON: 24-1.0103456-rdf.json
- Turtle: 24-1.0103456-turtle.txt
- N-Triples: 24-1.0103456-rdf-ntriples.txt
- Original Record: 24-1.0103456-source.json
- Full Text
- 24-1.0103456-fulltext.txt
- Citation
- 24-1.0103456.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0103456/manifest