Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Measuring functional connectivity of the brain
(USC Thesis Other)
Measuring functional connectivity of the brain
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
MEASURINGFUNCTIONALCONNECTIVITYOFTHEBRAIN by SergulAydore ADissertationPresentedtothe FACULTYOFTHEUSCGRADUATESCHOOL UNIVERSITYOFSOUTHERNCALIFORNIA InPartialFulfillmentofthe RequirementsfortheDegree DOCTOROFPHILOSOPHY (ELECTRICALENGINEERING) December2014 DoctoralCommittee: ProfessorRichardM.Leahy,Chair ProfessorJustinP.Haldar ProfessorKrishnaS.Nayak ProfessorVasilisZ.Marmarelis B.KeithJenkins Copyright 2014 SergulAydore Abstract The rich temporal content of measurements of electromagnetic activity, including elec- troencephalography (EEG) and magnetoencephalography (MEG), allow researchers to study dynamic functional networks in the human brain. However, it is difficult to turn this data into meaningful conclusions about those brain networks. In this dissertation, we describe the theoretical relationships between different interaction measures, fol- lowed by development of novel measures to address classical nuisance of cross-talk in brainelectrophysiologicalrecordings. Coherenceand phase lockingvalue(PLV)are widelyusedmeasures that can reveal interactions between electrophysiological signals within a frequency range of interest. We investigate the statistical properties of the PLV by describing two distributions that are widely used to a priori model phase interactions. The first of these is the von Mises distribution, for which the standard sample PLV is a maximum likelihood estimator. Thesecondistherelativephasedistributionderivedfrombivariatecircularlysymmetric complexGaussiandata. WederiveanexplicitexpressionforthePLVforthisdistribution and show that it is a function of the coherence between the two signals. We then com- pare results via local field potential data from a visually-cued motor study in macaque for the two different PLV estimators and conclude that, for this data, the sample PLV provides equivalent information to the coherence of the two complex time series. This ii resultreducestheanalysisoftime-lockedactivitybetweensignalstothecomputationof coherenceratherthancoherenceandPLV. Since the PLV is a bivariate measure (that is, it is computed pairwise between signals) it cannot differentiate between direct and indirect connections in a multidi- mensional network. A non-parametric partial phase synchronization index attempted to resolve this problem by extending sample PLV to the multivariate case using the same mechanism relating correlation to partial correlation. Here we derive an analyti- cal expression for partial PLV for a multivariate circular complex Gaussian model and show that partial PLV can be computed from partial coherence. We demonstrate our method in simulations with Roessler oscillators and experimental data of multichannel local field potentials from a macaque monkey. We show that the multivariate circular complex Gaussian model suggests similar synchronization networks. In addition, the circularcomplexGaussianmodelhasalowervarianceintheestimationofpartialPLV. Interpretation of functional connectivity from EEG/MEG data is challenging due to cross-talk problem between signals of interest. For example, coherence may yield spu- riously large values leading false positive connections. Approaches such as imaginary coherence,phaselagindex,laggedcoherenceandorthogonalcoherencehavebeenpro- posed to overcome this problem. The common assumption of these measures is that time-lagged interactions are more robust to cross-talk than instantaneous interactions. However,noneofthesemeasuresaccountfortheremainingnodesinamultivariatenet- work. While partial coherence quantifies the direct relationship between signals after excluding the linear effect from the remaining signals, it is still significantly affected by cross-talk. Here we combine the cross-talk robustness of lagged coherence with the multivariate framework of partial coherence to form a new measure called partial lagged coherence. Briefly, partial lagged coherence regresses our signals of interest iii onto the remaining signals so that only the instantaneous contributions from other sig- nals are removed before computing lagged coherence. Our findings on realistic sim- ulations of MEG data indicate a better performance of partial lagged coherence than other approaches in distinguishing direct from indirect connections in the presence of cross-talk. iv tomy mother and all the other courageouswomen in Turkey v Acknowledgements Iamgratefultohavehadinfinitesupportfromallthebeautifulpeoplearoundmeduring my doctoral studies over the years. This dissertation could not have been completed withoutthem. Iwouldliketoexpressmydeepestgratitudetothosewonderfulpeople. IamveryfortunateforbeingmentoredbyadistinguishedscientistProfessorRichard M.Leahyduringmydoctoralstudies. Heguidedmetoformulateimportantproblemsin thefieldwhileencouragingmetobeindependentovertime;hespentenormousamounts of time with me discussing not only science but also how to write research papers and howtopresentitinfrontofanaudience;healwaysencouragedmetoattendconferences and workshops which enhanced my training. It is hard to imagine a better mentor than him for a Ph.D. study. I would like to gratefully and sincerely thank Professor Richard Leahyforhisguidance,patience,encouragementandunderstanding. IwouldliketothankProfessorJustinP.Haldarforagreeingtoserveonthecommit- teeandhisintellectualcontributionstothedevelopmentofmyresearchduringourgroup meetings. I have incurred many debts to the rest of my committe members: Professor Krishna S. Nayak, Professor Vasilis Z. Marmarelis and B. Keith Jenkins for their par- ticipation, academic support and patience. I truly appreciate all of their time, guidance andsuggestionstowardimprovingthepreliminaryversionofthisdissertation. A very special thanks is due to my dear friend and colleague Dr. Syed Ashrafulla whomIconsideredasthebestfriendinLosAngelesforbeinganawesomefriendanda vi greatcollaborator. Hisexistencemademygraduateschoollifemorejoyfulandintellec- tually more productive. This dissertation would not have been possible without Syed’s supportandencouragement. I also want to give thanks to Professor Anand A. Joshi whom was the first person I talkedtowhenIneededalisteningearandadissentingopinion. Iparticularlyappreciate hisadviceonbuildingacareerinacademia. IhavebeenveryluckytohavehadgreatcollaboratorsasDr. JohnC.Mosher,Profes- sorSylvainBailletandFranc ¸oisTadelwhomademeawareofthepracticalimplications ofmymostlytheoreticalwork. I thank Professor Robert A. Scholtz who has been very kind to extend his help at variousphasesofmyresearchespeciallyinstatisticalanalysisofcomplexrandomvari- ables. The papers and the books he recommended and his lecture notes have been a greathelptowarddevelopmentofthechapters2and3. I was very privileged to be a member of the most friendly research group. I do, hereby,acknowledgeallthecurrentandpastmembersofthisawesomegroup: Dr. David Wheland, Dr. Dimitrios Pantazis, Dr. Joyita Dutta, Dr. Sangeetha Somaya- jula, Pro- fessor Juan Luis Poletti Soto, Dr. Yu Teng Chang, Chitresh Bhushan, Yanguang Lin, WentaoZhu,DivyaVaradarajanandDaeunKim. Thankyouallformakingtheresearch environmentsuchawarmandfunplace. I want to thank my former advisors and co-advisors Professor Kıvanc ¸ Mıhc ¸ak, Pro- fessor Ata Akın and Professor Yasemin Kahya from Bo˘ gazic ¸i University who encour- agedandsupportedmetopursueanacademiccareeraftermymasterdegree. I owe a great dept of gratitute to my life-long friends Filiz Ates ¸, Nemin Topalo˘ glu, Helin Duta˘ gacı, Berc ¸em Duta˘ gacı who have always been there despite the hundreds of miles between us. I am eternally grateful to my dear friends Deniz C ¸akırer, Ekin Pehli- van, Yaman ¨ Ozakın for their caring and support especially during the most challenging vii phase of my doctoral studies. Thank you all for bringing courage, wisdom and joy to mylife. Duringmydoctoralstudies,Ifoundaperfectwaytoescapefromaccumulatedstress ofmydailylifethroughTango-themostsensualandbeautifulformofdance. Iwantto thankmygreatTangoinstructorsandfriendsMitraMartinandStefanFabryforhelping me discover the harmony between myself and the blue planet. I also want to thank all thedancersatOxygenTangofortheirwarmembraces. Above all, I owe a lot to my mother Esen Ayd¨ ore and my father Serdar Ayd¨ ore for their love, support and sacrifices in every stage of my life. They have taught me how to be persistent, independent and hard working. I have managed to deal with all the challenges I have faced during my doctoral studies through their unwavering belief in me. I saved the last word of the acknowledgements for my brother Onur Ayd¨ ore who is the most important person in my life. Thank you for giving me many happy and beautifulmemoriesthroughoutmylifeandbeingmybiggestsupporter. viii Contents Abstract ii Dedication v Acknowledgements vi ListofFigures xii ListofTables xvii 1 Introduction 1 1.1 TheHumanBrain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 History . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.2 Anatomy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2 FunctionalBrainImaging . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.2.1 PositronEmissionTomography(PET) . . . . . . . . . . . . . . 11 1.2.2 Functionalmagneticresonanceimaging(fMRI) . . . . . . . . . 13 1.2.3 Electrocorticography,LocalFieldPotentialsandSingle-neuron ActionPotentialRecordings . . . . . . . . . . . . . . . . . . . 14 1.2.4 EEGandMEG . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.2.5 ComparisonofImagingTechniques . . . . . . . . . . . . . . . 20 1.3 ConnectivityintheBrain . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.3.1 StructuralConnectivity . . . . . . . . . . . . . . . . . . . . . . 23 1.3.2 FunctionalConnectivity . . . . . . . . . . . . . . . . . . . . . 24 1.4 BrainDisorders . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.4.1 Schizophrenia . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.4.2 MajorDepression . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.4.3 Alzheimer’sdisease . . . . . . . . . . . . . . . . . . . . . . . 29 1.4.4 Epilepsy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 1.5 OrganizationoftheDissertation . . . . . . . . . . . . . . . . . . . . . 30 ix 2 PhaseSynchronization 31 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.2 MeasuresofPhaseSynchronization . . . . . . . . . . . . . . . . . . . 33 2.2.1 ThePhaseLockingValueandPhaseLagIndex . . . . . . . . . 33 2.2.2 PhaseLockingValueandthevonMisesDistribution . . . . . . 37 2.2.3 Phase Locking Value and Circularly Symmetric Gaussian Pro- cesses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.3 Results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 2.3.1 SimulationsBasedontheGaussiandistribution . . . . . . . . . 45 2.3.2 RoesslerOscillatorSimulations . . . . . . . . . . . . . . . . . 46 2.3.3 AnalysisofLFPdata . . . . . . . . . . . . . . . . . . . . . . . 49 2.4 DiscussionandConclusion . . . . . . . . . . . . . . . . . . . . . . . . 54 3 PartialPhaseSynchronization 60 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.2.1 Non-parametricPartialPLV . . . . . . . . . . . . . . . . . . . 61 3.2.2 Parametric partial PLV Using the Multivariate Circular Com- plexGaussianModel . . . . . . . . . . . . . . . . . . . . . . . 62 3.3 Results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.3.1 RoesslerOscillatorSimulations . . . . . . . . . . . . . . . . . 64 3.3.2 Localfieldpotentialsfromamacaquemonkey . . . . . . . . . 67 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4 ConnectivityinthePresenceofCross-talk 70 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.2 ProblemDefinition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.3.1 InteractionMeasuresinPresenceofaSingleCross-talk . . . . . 75 4.3.2 InteractionmeasuresinPresenceofmultipleCross-talk . . . . . 80 4.4 Results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.4.1 SingleCross-talk . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.4.2 MultipleCross-talk . . . . . . . . . . . . . . . . . . . . . . . . 85 4.5 DiscussionandConclusion . . . . . . . . . . . . . . . . . . . . . . . . 93 5 ConclusionandFutureDirections 97 5.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.2 FutureDirections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2.1 PartialPhaseSynchronization . . . . . . . . . . . . . . . . . . 99 5.2.2 ConnectivityinthePresenceofCross-talk . . . . . . . . . . . . 100 5.2.3 ClinicalApplications . . . . . . . . . . . . . . . . . . . . . . . 100 x ReferenceList 102 A PhaseSynchronization 119 A.1 BiasofSquaredPLV . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 A.2 UnbiasedsquarePLVandPPC . . . . . . . . . . . . . . . . . . . . . . 119 A.3 Conditional Distribution of Relative Phase for the Bivariate Circularly SymmetricComplexGaussian . . . . . . . . . . . . . . . . . . . . . . 120 A.4 ThePLVfortheBivariateCircularlySymmetricComplexGaussian . . 121 B PartialPhaseSynchronization 126 B.1 Conditionalcovariancematrixofxconditionedony . . . . . . . . . . 126 C ConnectivityinPresenceofLinearMixing 129 C.1 Linear Mixing Invariant Property of Lagged Coherence in presence of twosources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 C.2 Proof of the relationship between Orthogonal Coherence and the Stan- dardCoherence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 xi ListofFigures 1.1 SketchesofthehumanskullandbrainbyLeonardodaVinci. (a)Draw- ingofahumanskull(WindsorCastle,RL19058recto,c.1489)(b)Two drawings of skulls (Windsor Castle, RL 19057 recto, c. 1489) . The rectangle drawn around the lower skull formed part of Leonardos stud- ies of the proportions of the human body. (c) Central nervous system and cranial nerves (Windsor Castle, RL 12603 recto, c. 1490). The maindrawingshowsthelayerscoveringthebrain,comparedtothelay- ersofanonion. Themaindrawingandthesketchbelowshowthethree cerebralventricles,followingthetraditionalmedievalframeworkofIbn Sinaandothers. (d)Twodrawingsofthehumanhead(detailofadraw- ingfromtheSchlossmuseum,Weimardatedfromc.15061507). Figures and captions are from [Pev02]. A translation of the text is provided at http://pevsnerlab.kennedykrieger.org/leonardo.htm. . . . . . . . . . . . 2 1.2 Three main parts of the human brain: the forebrain, the midbrain and thehindbrain. Imagecredit: http://www.studyblue.com . . . . . . . . . 5 1.3 Fourlobesofthecerebralcortex. Imagecredit: http://www.deryckthake.com 7 1.4 The structure of the human cerebral cortex, including the association cortices. (A)Asummaryofthecellularcompositionofthesixlayersof the cerebral cortex. (B) Based on variations in the thickness, cell den- sity, and other histological features of the six neocortical laminae, the humanbraincanbedividedintoabout50cytoarchitectonicareas,inthis casethoserecognizedbytheneuroanatomistKorbinianBrodmanninhis seminalmonographin1909. Fromhttp://www.ncbi.nlm.nih.gov/books/NBK10870/ 8 1.5 Diagramofaneuronandsynaptictransmission. Imagecredit: http://www.brainfacts.com 9 1.6 SchemaofaPETacquisitionprocess. Imagecredit: http://www.wikipedia.com 12 1.7 Schema of the change in oxyHb and deoxyHb during neuronal activity. Bloodoxygenationincreaseswithanincreasedactivitybutonlyamod- est amount of oxygen is consumed which results in a decrease in the concentrationofdeoxyHb. From[HSM04] . . . . . . . . . . . . . . . 14 1.8 Implantationofintracranialelectrodes. Imagecourtesy: ClevelandEpilepsy Clinic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 xii 1.9 ComparisonoftheSpatialDomainsofthethreeinvasiveElectricalRecord- ing Modalities and noninvasive EEG. A typical noninvasive EEG elec- trode(blackcircle)locatedapproximately2cm2abovethecortexmea- sure average gyral neuronal activity across a 3 cm spatial extent (filled black cortical layers). ECoG electrodes lie just above the cortex and average neural activity over a smaller 0.5 cm range. Both local field potentials and single-unit action potentials are recorded fromwithin the brainparenchymaandsampleevensmallerareasofneuraltissue,yield- inghigherspatialresolutions. Figurefrom[SCWM06]. . . . . . . . . . 16 1.10 Depiction of the basics of EEG and MEG. (a) An ionic current gen- erated by large neuronal cells. The primary current is represented by a black large arrow while the dark lines represent secondary currents. Magneticfieldswhicharegeneratedbyprimaryandsecondarycurrents are represented by dashed lines. (b) The neuronal populations generate mass effect of currents which behave as a current dipole (red arrow). EEG detects the secondary currents (circular yellow lines) induced by primarygenerator. MEGcapturesthemagneticfieldswhichtravelmore easilywithintissues. From[Bai12]. . . . . . . . . . . . . . . . . . . . 18 1.11 (a) An elastic EEG cap with 60 electrodes. From [Bai12]. (b) EEG electoresrelativetothehumanheadandthefoldedcortex. From[Nun95]. 18 1.12 A typical MEG equipment. MEG measurements are obtained inside a magnetically shielded room (left). Liquid helium (upper right) stored in a Dewar (left and upper right) keeps MEG sensors which use low- electronic temperature cool. Visualization of time-evolving magnetic fieldtopographies(lowerright). From[BML01]. . . . . . . . . . . . . 19 1.13 Structural and Functional Brain Networks. First, the network nodes are defined for both connectivity types. Second, measures of interaction bothforstructuralandfunctionalconnectivityarecomputed. Thiscould bethecorrelationmetricsbasedontemporalinteractionssuchascoher- ence, Granger causality for functional networks while the connection probabilitybetweentworegionsofanindividualdiffusiontensorimag- ingdataset,ortheinter-regionalcorrelationsincorticalthicknessorvol- ume MRI measurements estimated in groups of subjects for structural networks. Finally, connectivity networks are obtained after statistical thresholding. From[BS09]. . . . . . . . . . . . . . . . . . . . . . . . . 22 2.1 Bias and variance of ^ PLV 2 sample (red) and ^ PLV 2 ub sample (blue) as a function of N, the number of samples used to compute the estimators: (a) mean value vs N for four different true values of PLV; (b) variance vs. N for the same true values of PLV. Samples were drawn indepen- dently from the von Mises distribution for four different concentration parametervalues. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 xiii 2.2 Probabilitydensityfunctions: (a)vonMisesdistributionwithmean = 0 for a range of concentration values, ; (b) relative phase distribution for bivariate circularly symmetric Gaussian models for different values ofcross-correlationmagnitudejR 12 j = 12 p 11 22 withphase 12 = 0. . . 38 2.3 Plots of sample (a) PLV and (b) PLI as a function of concentration parameter and mean for samples drawn from the von Mises dis- tribution. Note that while PLV is independent of , PLI depends on bothand. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.4 Plot of the monotonic relationship between PLV and the magnitude of cross-correlation,jR 12 j = 12 p 11 22 ,forthebivariatecircularlysymmetric Gaussianmodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.5 Plots of sample (a) PLV and (b) PLI as a function of cross-correlation magnitude jR 12 j = 12 p 11 22 and phase 12 for the bivariate circularly symmetric Gaussian distribution. Similarly to the von Mises distribu- tion, PLV is independent of phase while PLI depends on both cross- correlation magnitude and phase. Samples were drawn from relative phasedistributioninequ. (2.22). . . . . . . . . . . . . . . . . . . . . . 44 2.6 Plot of bias and variance for the sample PLV, root of unbiased sample PLV, and PLV computed from sample coherence. (a) plot of the mean over 1,000 Monte Carlo trials as a function of number of samples used toestimatetheparameterfortwodifferentvaluesofcoherence;(b)cor- responding plot of variance for each of the three estimators for the two different coherence values. Samples were drawn independently from bivariateGaussianprocesses. . . . . . . . . . . . . . . . . . . . . . . . 45 2.7 Roessleroscillatorsimulations(a)ComparisonofROCcurvesof ^ PLV circgauss (blue)and ^ PLV sample (red)whenϵ = 0:15forcoupledoscillators,stan- dard deviation = 1:5 and number of samples L = 5000. (b) Area under ROC curve as a function of coupling parameter ϵ when = 1:5 and L = 5000. (c) Area under ROC curve as a function of number of samples L when = 1:5 and ϵ = 0:15. (d) Area under ROC curve as a function of the standard deviation of the noise whenL = 5000 and ϵ = 0:15. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 2.8 Locationsof15electrodepairsintherighthemisphere(reproducedfrom [LDB01]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 2.9 Four visual cues represented to the monkey in the experiment (a) right slantedline(b)leftslantedline(c)rightslanteddiamond(d)leftslanted diamond. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 xiv 2.10 Goodnessoffitbetweenempiricaldistributionofrelativephaseextracted fromLFPdataandparametricvonMisesandcircularlysymmetricGaus- sian distributions. The concentration parameter (von Mises) and cross- correlationparameter(Gaussian)weredirectlyestimatedfromthesam- pledata. (a-c)pairsat 120msec(d-h)pairsat260msec . . . . . . . . . 58 2.11 Scatterplotof ^ PLV vonmises versus ^ PLV circgauss computedfrommacaque LFPdata. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.12 Phase synchronization networks constructed from the results in Table 2.2. Green: significantPLVsoccurwhendiamondispresented. Magenta: significant PLVs occur when line is presented. Blue: significant PLVs occur during go condition. Red: significant PLVs occur during no-go condition. The thickness of the edges represents the number of experi- mentssessionsinwhichsignificantsynchronizationoccurs. . . . . . . . 59 3.1 (a) Estimated coupling between the 2nd and 3rd Roessler oscillators. For each boxplot, the central mark is the median, the edges of the box are the 25-th and 75-th percentiles, and the whiskers extend to the most extreme data points not considered outliers. The true coupling between is zero. (b) Q-Q plot of ranked Mahalanobis distance r 2 i versus the expectedcorrespondingvalueof 2 3 fortheRoessleroscillators. . . . . . 65 3.2 ROCanalysisoftheRoessleroscillators(a)ROCcurvesforϵ = 0:2(b) AreaunderROCcurvesfordifferentvaluesofϵ(legendas(a)). . . . . . 66 3.3 Q-QplotofrankedMahalanobisdistanceversusthecorrespondingquan- tileof 2 15 formacaqueLFPdata. . . . . . . . . . . . . . . . . . . . . . 68 3.4 (a-c) Bootstrap mean of the PLV measures. (d) Difference of bootstrap standarddeviationofnon-parametricminusparametricpartialPLV. . . 68 4.1 Differentinterferenceandlinearmixingschemes. Thesourcess 1 ands 2 aremeasuredinchannelsx 1 andx 2 . (a)Idealcase: nolinearmixingno interference. (b) linear mixing without interference. (c) linear mixing withinterference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.2 Interaction measures as a function of cross-talk coefficients in presence ofasinglecross-talk. (a)Coherence(b)Imaginarycoherence(c)Phase lag index (d) Lagged coherence. Except from lagged coherence (d), all other interaction measures have high dependence on the amount of cross-talk. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4.3 Eightextendedsourcesareselectedonthecorticalsurfacebasedonthe locationsofimplantedelectrodesinamacaquemonkeyexperiment. . . 87 xv 4.4 Average values of measures over 100 Monte Carlo simulations where only the sources 2 and 7 are connected. Since the measures are sym- metric only the upper halves of the connectivity matrices are shown. (a)Coherence(b)ImaginaryCoherence(c)PartialCoherence(d)Phase LagIndex(e)LaggedCoherence(f)PartialLaggedCoherence . . . . . 89 4.5 NormalizedHistogramsofconnectedandunconnectededges . . . . . . 90 4.6 ComparisonofROCcurvesofdifferentinteractionmeasuresinpresence ofasingleconnectioninamultivariatenetwork. . . . . . . . . . . . . . 92 4.7 Truenetworksbasedonarandomselection . . . . . . . . . . . . . . . 93 4.8 FractionalareaunderROCasafunctionofnumberofconnections . . . 94 xvi ListofTables 2.1 Thep-valuesofChi-squarestatisticforgoodnessoffitofthevonMises and circularly symmetric Gaussian phase models. In only on case are we able to reject the null hypothesis: electrode pair 3 and 5 for the von Misesmodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 2.2 Comparison of the number of experiments showing significant interac- tionsbetweenelectrodepairsforeachoffourconditions(Diamond,line, go, no-go) computed using sample PLV and Gaussian PLV for early and late response periods. Significance was assessed using FDR = 0.01 basedonpermutationsoftrialindices. Onlypairsforwhichaminimum of5experimentsshowsignificantinteractionsineithercueortaskbased comparisonareshown. . . . . . . . . . . . . . . . . . . . . . . . . . . 53 xvii Chapter1 Introduction 1.1 TheHumanBrain The brain - the most complex organ in the human body- has fascinated scientists and philosophers for centuries. This organ is responsible for our thoughts, actions, emo- tions, perceptions and memories which make each of us unique. Thanks to the recent advanced technologies and developments of brain research, we have started to unlock the complexity of the brain. This section is a basic introduction to the brain and recent brainimagingtechnologies. 1.1.1 History Thehistoryofhumanbrainstudiesisasoldasthehumanhistoryitself[Ade87]. Archae- logicalfindingssuggestthepracticeofprehistorichumans’primitivebrainsurgeriesfor treatmentsinSouthAmerica. However,hearthadbeenthoughttobethecentralorganof ourbodyforcenturies untiltheearly GreekphysicianAlcmaeonsuggested brainasthe central organ of our thoughts and sensations (around 450 B.C.). Based on this theory, Hippocrates believed that the brain is the seat for intelligence. The Greek philosopher Aristotles, on the other hand, stated that heart is the seat of intelligence and the brain wasacoolingmechanismforthebloodbuthesuccessfullystatedadistinctionbetween 1 Figure 1.1: Sketches of the human skull and brain by Leonardo da Vinci. (a) Draw- ing of a human skull (Windsor Castle, RL 19058 recto, c.1489) (b) Two drawings of skulls (Windsor Castle, RL 19057 recto, c. 1489) . The rectangle drawn around the lower skull formed part of Leonardos studies of the proportions of the human body. (c) Central nervous system and cranial nerves (Windsor Castle, RL 12603 recto, c. 1490). The main drawing shows the layers covering the brain, compared to the layers of an onion. The main drawing and the sketch below show the three cerebral ventricles, fol- lowing the traditional medieval framework of Ibn Sina and others. (d) Two drawings of the human head (detail of a drawing from the Schlossmuseum, Weimar dated from c.15061507). Figures and captions are from [Pev02]. A translation of the text is pro- videdathttp://pevsnerlab.kennedykrieger.org/leonardo.htm 2 short term memory and long term memory (around 335 B.C.). Two Alexandrian biolo- gists Herophilus and Erasistratus discovered the human nervous system by dissecting a humanbody(around330B.C.). Galen,aninfluentialphysicianduringRomanEmpire, proposed that our emotions and functioning of our body are determined by the fluids andhumorsinthebrain[Ade87]. Unfortunately, brain studies in Europe ceased during the middle ages due to a ban on the study of anatomy by the church while classical medicine continued to develop in middle Eastern civilizations. Hence, the famous writings of medical scholars arrived EuropethroughSyrian,Hebrew,PersianorArabictranslationsintoLatin[Ade87]. The idea of naturalism in art during the Renaissance induced an interest in human anatomy. Dissection of human cadavers were encouraged in public and attracted large audiences. Andreas Vesalius, a Renaissance anatomist, reported many anatomical fea- turesofthebrainduringhisdissectionsofhumancadavers(around1543). Hepublished his findings in the first known neuroscience textbook: “On the fabric of human brain”. Heopposedthetheorythatthehigherfunctionsofthebrainaresituatedintheventricles after his discovery of the existence of the same ventricles in animals [VL92]. His work was extended by Leonardo Da Vinci and depicted in a detailed sketches of the human brain[Pev02]. Fig. 1.1illustratestheintersectionofhisextraordinaryartandscience. Rene Descartes (1596-1650), a French mathematician and a philosopher, proposed a distinction between immaterial “mind” and organ of the brain based on his theory of dualism. He argues that a “soul” is contained in the mind not in the brain [Ade87]. He furthersuggestedthatthesetwointeractatthepinealgland. ThomasWillis(1621-1675), anatomist and physician at Oxford, illustrated the most detailed description of brain’s anatomy to that date in his Cerebri Anatome. He stated that the thoughts and voluntary movementwerecontrolledbythecerebrumwhichconstitutes70%ofthebrain. Reflex movement, on the other hand, was controlled by the cerebellum. Luigi Galgavani, an 3 Italianphysiologist,pioneeredthestudiesofelectrophysiologythroughhisexperiments on a frog. He interpreted his results as the flowing property of electricity from muscles tothenervewhichlaterdisprovedbyAlessandroVolta(1745-1827). In 1848, a railroad worker Phineas Gage was involved in a bizarre accident where an iron bar shot through his head destroying his frontal lobe [Mac02]. Despite his sur- prisingsurvival,hischaracterhadprofoundlychangedfromaquietmantoacombative man. This incident was a milestone in neuroscience suggesting the key role of frontal lobeinapersonality. WhathappenedtoGageplayedanindirectroleinthedevelopment ofaprocedurecalledlobotomywhichconsistsofremovingportionsofthefrontallobe. PaulBroca(1824-1880),aphysician,surgeonandresearcher,discoveredthespeech center of the brain in the frontal lobe which has been named after him as Broca’s area [GA06]. Inparallelwiththisresearch,CarlWernicke(1848-1905)publishedhisdiscov- eryofthespecializedareainthebrainoflanguagecomprehensionandproductioninhis essay Der Aphasische Symptomencomplex. This part of the brain is called Wernicke’s area[HM77]. Camillo Golgi (1843-1926), an Italian scientist, developed a staining procedure to reveal the structures of single neurons [BK68]. Santiago Ramon y Cajal (1852-1934), a Spanish neuroscientist, applied his method to formulate the hypothesis that the basic functional unit of the brain is the neuron [She91]. Golgi and Cajal shared the Nobel PrizeinPhysiologyorMedicinein1906fortheirobservationsonneurons[Gra07]. Advances in the functional brain imaging technologies marked the 20-th century. Newbrainmappingmodalitiesallowedresearcherstomonitortheactivationinthebrain. Wewillintroducethesemodalitiesinsection1.2. 4 Figure 1.2: Three main parts of the human brain: the forebrain, the midbrain and the hindbrain. Imagecredit: http://www.studyblue.com 1.1.2 Anatomy Thebrainisoneofthelargestandmostimportantorgansofthehumanbody. Thisthree pound organ consists of three main parts: the forebrain, the midbrain and the hindbrain asillustratedinFig. 1.2. Themidbrainandthehindbrainarealsoknownascentralcore and it is differentthan the forebrain. The central core is relativelysimple and wasquite static during the evolutionary process of the human brain. It is largely responsible for unconscious activities. The forebrain, on the other hand, is highly complex and is the sourceforintellectualactivities[DCB91]. The hindbrain is also called the brainstem. It connects the rest of the brain to the spinalcordandincludesthemedullaandthecerebellum. Themedullacontrolsthevital functions such as respirations, digestion and blood circulation. Cerebellum controls movement,balanceandcoordination. Itisalsoinvolvedinlearnedactivities[DCB91]. 5 Themidbrainislinkedtohumanemotionandlongtermmemoryviathelymbicsystem it contains. The lymbic system is located above the brain stem and under the cerebral cortex and very similar in all mammals. The features included in the lymbic system are: the hypothalamus, the amygdala and the hippocampus. Physical changes in heart rateandrespirationaccompaniedbystrongemotionsareregulatedinthehypothalamus. The amygdala is associated with aggression and emotion. The hippocampus is the key featuretolongtermmemory[DCB91]. Theforebrain isthe largestand mostdevelopedpart ofthe humanbrain. It contains cerebral cortex that covers the brainstem and the midbrain. The cerebral cortex gives humansabilitytosolvecomplexproblems. Itisthepartthatdistinguisheshumansfrom other animals [DCB91]. This dissertation focuses on the functions of cerebral cortex. Therefore,itwillbeexplainedinamoredetailedwayinthefollowingsection. TheCerebralCortex The cerebral cortex, a large folded surface, is the outer mantle of cells of the human brain which is believed to play an important role for the brain to process higher-level functions [TE98, Hor03]. This is the part we see when we look at the human brain. The cerebral cortex has evolved into two cerebral hemispheres where usually one is dominant. Forinstance,leftcerebralhemisphereisdominantinrighthandedpeopleand vice versa. The folded structure creates structures called gyri (ridges on the cerebral cortex)andsulci(groovesinthecerebralcortex). Four lobes with different specialized functions make up the cerebral cortex. These lobesare: thefrontallobe,theparietallobe,thetemporallobeandtheoccipitallobe(see Fig. 1.3) [JP70]. The frontal lobe is located in front part of the brain and is responsible for reasoning, decision making, voluntary muscle movements, processing speech and high level cognition. The parietal lobe is located in the middle part of the brain and 6 is responsible for processing sensory information such as taste, temperature and smell. The temporal lobe is located in the bottom part of the brain and plays a key role in processing auditory stimuli. The occipital region is located in the back of the brain and heavilyassociatedwithinterpretingvisualinformation. Figure 1.3: Four lobes of the cerebral cortex. Image credit: http://www.deryckthake.com The six layered structure represents the great majority of the cerebral cortex where eachlayercomprisesdifferentcellpopulationsbasedontheirsizes,shapesandfunctions [cA]. Fig. 1.4illustratesthelayeredstructureandhistologicallydefined50subdivisions ofthecerebralcortex. These six layers are numbered with Roman numerals from outside to inside. Layer I is the molecular layer, which contains few scattered neurons; layer II, the external granular layer, contains small pyramidal neurons and numerous stellate neurons; layer 7 Figure 1.4: The structure of the human cerebral cortex, including the association cor- tices. (A)Asummaryofthecellularcompositionofthesixlayersofthecerebralcortex. (B) Based on variations in the thickness, cell density, and other histological features of the six neocortical laminae, the human brain can be divided into about 50 cytoarchitec- tonicareas,inthiscasethoserecognizedbytheneuroanatomistKorbinianBrodmannin hisseminalmonographin1909. Fromhttp://www.ncbi.nlm.nih.gov/books/NBK10870/ III ,the external pyramidal layer, contains predominantly small and medium-size pyra- midal neurons; layer IV, the internal granular layer, contains different types of stellate and pyramidal neurons; layer V, the internal pyramidal layer, contains large pyramidal neurons ; and layer VI, the multiform, or fusiform layer, contains few large pyramidal neurons and many small spindle-like pyramidal and multiform neurons. Each cortical layercontainsdifferentneuronalshapes,sizesanddensityaswellasdifferentorganiza- tionsofnervefibers. 8 TheNeuron The neuron (nerve cell) is the basic functional unit of the nervous system. An aver- age human brain contains about 100 bilion neurons. Neurons have ability to transmit informationtootherneuronsbyelectro-mechanicalsignalling[cA]. A typical neuron possesses a cell body, dendrites and a single axon as illustrated in Fig. 1.5. The cell body keeps the neuron alive and functioning. The dendrites extend fromthecellbodyandreceiveinputsfromotherneurons. Theaxonisasinglebranching structure extending from the cell body and ending in a number of synapsis (a junction between cells) where the signal can be delivered. It can be quite long compared to the other parts of the neurons. Most of the axons are covered with a myelin sheath to help acceleratingthesignaltransmissionalongtheaxon[cA]. Figure 1.5: Diagram of a neuron and synaptic transmission. Image credit: http://www.brainfacts.com 9 The communication between neurons is utilized via electrical signals and chemical messengers An information from other neurons or sensory receptors is received by a neuron through the dendrites. This piece of information travels to the cell body and then on to the axon in the form of an electrical signal which is called action potential. When the electric pulse reaches the synapses, it causes the release of chemical messen- gers called neurotransmitters, hence changes the electrical activity in the other neuron. Neurotransmitters are required to carry these electric pulses through the synaptic gap betweenneurons. Whenalargenumberofneuronsareactivatedinthisway,weareable tomeasurethisactiviationthroughelectromagneticdevices[cA]. Electrophysiology Action potentials also known as ”nerve impulses” or ”spikes” refer to electric signals generated by neurons. Some of the neurotransmitters released by a presynaptic cell are boundedtoreceptorsonthepostsynapticcell. Thisbindingopensvariousionchannels, allowing ions to enter or leave the cell. These ions are sodium (Na+), potassium (K+) and chloride (Cl-). The movement of electrically charged ions through ion channels accross the neuronal membrane alters the membrane potential. The synapse is excita- tory if the membrane is depolarized and inhibitory if the membrane is hyperpolarized. An excitatory postsynaptic current is generated at the excitatory synapses when the ion channel allows sodium into the cell. This current generates an excitatory postsynaptic potential (EPSP). ESPSs resulting from a single synapse is not sufficient to bring the neuron’s potential to its firing threshold. However, combined EPSPs from hundreds of synapses can cause a large fluctuations in membrane potential. Action potential is initiated when the membrane depolarization reaches a certain point called a threshold [Glo85]. It is conducted along the axon to cause neurotransmitter to be released at a 10 synapse on another cell [Luc17]. This mechanism plays a central role in communica- tionbetweenneurons. 1.2 FunctionalBrainImaging The electrical and magnetic fields generated by neurons that were described in section 1.1.2canberecordedbymeansofelectrodesatashortdistancefromthesources(single unit, LFP), or from cortical surface (ECoG) or at longer distances (EEG, MEG). The increase in neuronal activity also results in an increase in cerebral blood flow (CBF) to the active region due to the increased metabolic demand for glucose and oxygen. These hemodynamic changes can be detected via functional imaging modalities such as positron emission tomography (PET) and functional magnetic resonance imaging (fMRI). Functionalneuroimagingtechniquesinvolvesfunctionalandstructuralidentification oftheactivityofthebrain. Wedonotuseallthesetecniquesinthisdissertationbutitis crucial to understand what these modalities offer and their strengths and limitations in order to interpret the current problems in neuroimaging. The purpose of this section is toprovideabriefintroductiontocommonlyusedneuroimagingtechniques[CFM + 10]. 1.2.1 PositronEmissionTomography(PET) Positron emission tomography (PET) is a nuclear medical imaging technique that can noninvasivelymonitorthedistributionandreal-timemovementofinjectedradioactively labeled metabolically active chemicals, called radiotracers. When the tracer makes its waytothebrain,thesensorsinthePETscannerdetectsthegammaraysemittedthrough positron-electron annihilation in different regions of the brain (see Fig. 1.6). Three 11 dimensional images of the distribution of the chemicals can be produced via comput- erized reconstruction from the emission data. There are various tracers available but fludeoxyglucose (FDG), an analogue of glucose, is the most commonly used. In this case, PET images of the brain reflect the activity of the brain in terms of oxygen and glucosemetabolism[BTVM05]. Figure 1.6: Schema of a PET acquisition process. Image credit: http://www.wikipedia.com ThebiggestchallengeofPETscanningistherapiddecayofradioactivitywhichprevents researchers from conducting long experiments. Another drawback of PET scanning is theexposureofsubjectstoradiationandhencethelimitationinthenumberofscansthat thisimplies. 12 1.2.2 Functionalmagneticresonanceimaging(fMRI) Functional magnetic resonance imaging (fMRI) is a noninvasive technique to measure brainactivityusingMRIphysics. Thistechniquereliesonthelinkbetweenneuralactiv- ityandcerebralbloodflowinthebrain. Anincreasedactivityinaparticularbrainregion demands more oxygen and hence more blood flow. The molecule called hemoglobin in red blood cells carries oxygen to neurons. fMRI utilizes the magnetic property of hemoglobin which is paramagnetic when deoxygenated (deoxyHb) and diamagnetic whenoxygenated(oxyHb). ThevesselsthatcarryoxyHbdonotcauseasignificantdis- tortion to the magnetic field while capillaries and veins that carry deoxyHb distorts the local magnetic field [PC36]. Inhomogenity of the magnetic difference in various parts ofthebrainisdetectedbywaterprotons. Thisaffectstherelaxationtimeofwaterwhich will modify MRI signals. The presence of deoxyHb shortens the relaxation time and leads to reduction in MRI signals. After a momentary decrease in the blood oxygena- tion, OxyHb increases and deoxyHb decreases with an increased neuronal activity (see Fig. 1.7). Therefore,relaxationtimebecomeslongerandhencetheMRIsignalintensity increases with respect to the baseline. Mapping the magnetic changes of deoxyHb this wayiscalledbloodoxygenationleveldependent(BOLD)imaging[HSM04]. ThespatialresolutionoffMRIisintheorderofafewmillimeters,whereasthetem- poralresolutionisaroundafewsecondduetotheslowvariationinthebloodoxyenation level[HSM04]. 13 Figure 1.7: Schema of the change in oxyHb and deoxyHb during neuronal activity. Blood oxygenation increases with an increased activity but only a modest amount of oxygenisconsumedwhichresultsinadecreaseintheconcentrationofdeoxyHb. From [HSM04] 1.2.3 Electrocorticography , Local Field Potentials and Single- neuronActionPotentialRecordings Theelectricalactivityincerebralcortexcanberecordedinvasivelybyplacingelectrodes on the exposed surface. Currently, the three primary recording modalities are electro- corticography (ECoG), local field potentials (LFPs), and single-neuron action potential recordings (single units). These techniques are classified by the location of the elec- trodesthatmeasureextracellularpotentialsgeneratedbyneuronsincorticallayer. ECoG uses subdural platinium-iridium or stainless steel electrodes to record elec- trical activity directly on cerebral cortex, thereby avoids the signal distortion caused by intermediate layers [BAK12] as shown in Fig. 1.8. The spatial resolution of ECoG is on the order of mm [GG01] in addition to the ability to record high frequency content inneuronalsignals[PGH + 03]. Ithasbeenperformedinclinicalsettings,particularlyin epilepsypatientstolocalizeseizureonsetzones[EMFO05]. 14 Figure1.8: Implantationofintracranialelectrodes. Imagecourtesy: ClevelandEpilepsy Clinic The low-pass filtered (usually with a cutoff< 250 Hz) electric potential occured in theextracellularspacearoundneuronsisalsocalledthelocalfieldpotential(LFP).This signal reflects the electrical events at deeper locations and can be recorded in the extra- cellular space in brain tissue, typically using micro-electrodes (metal, silicon or glass micropipettes). LFPs sample the activity from relatively smaller number of neurons compared to ECoG which uses larger electrodes implanted in the surface. Electrodes separated by 1 mm [DCS99] or by a few hundred microns [KNB + 09] can capture and distinguish the activity in different localized populations of neurons. The rich tempo- ral and spatial information content of LFPs opens up the possibility to investigate the intracellulardynamicsinthebrain. 15 Figure1.9: ComparisonoftheSpatialDomainsofthethreeinvasiveElectricalRecord- ing Modalities and noninvasive EEG. A typical noninvasive EEG electrode (black cir- cle) located approximately 2 cm2 above the cortex measure average gyral neuronal activity across a 3 cm spatial extent (filled black cortical layers). ECoG electrodes lie just above the cortex and average neural activity over a smaller 0.5 cm range. Both localfieldpotentialsandsingle-unitactionpotentialsarerecordedfromwithinthebrain parenchyma and sample even smaller areas of neural tissue, yielding higher spatial res- olutions. Figurefrom[SCWM06]. Single unit recordings sample action potentials from single, firing neurons with micro- electrodes which typically have 10 50 m tips, shanks that are 1 5 mm in length, and one or more conductors surrounded by insulation [LPA + 01]. Since the neurons are the basic functional units in the brain, single unit recordings provide the most pre- cise measures of activity in the brain. Single unit recordings have been used to study behaviors and functions in the brains of animals and humans [BPS + 99]. Recently, sin- gle unit recordings have been performed in brain machine interface studies which are used to control the movement of a device based on the decoded signals from the brain [SCWM06]. 16 ThecomparisonofthespatialdomainsofthesetechniquesandscalpEEGisshown in Fig. 1.9. In this dissertation, we restricted our attention to LFPs among the invasive techniques. 1.2.4 EEGandMEG Noninvasive measurement of the electric potential generated by a single neuron is not achievable. However, when thousands or millions of neurons fire synchronously, the summationofelectricandmagneticfieldscanbemeasuredoutsidethehead. MEGand EEG measure, respectively, magnetic and electric fields outside the head produced by neural cell assemblies. It is believed that the pyramidal neurons are the main sources of EEG and MEG because of their coherent distribution of their large dendritic trunks locallyorientedinparallelandpointingperpendicularlytotheconvolutedcorticalsheet ofgraymatter[HHI + 93,Sil95]. Thesetwoimagingmodalitiesofferauniquetemporaltimeresolutionbelow100ms. This is essential for resolving rapid change of neurophysiological process. Thus, we can study the dynamics of neural processes that occur at typical time scales on the order of of tens of milliseconds [Nun06]. Therefore, MEG and EEG have become two complimentary techniques for noninvasive assessment of the electrical activity in the brain. Figure1.10illustratesmeasurementoftheactivitygeneratedbylargeneuronalcells. Electric potential differences between pairs of electrodes on scalp are recorded in EEG as shown in Fig. 1.11(a). The positioning of the electrodes relative to the layers of the head are depicted in Fig. 1.11(b). In contrast, magnetic fields produced by neuronal currents is very weak compared to the EEG scalp voltages. Thanks to the invention 17 (a) (b) Figure 1.10: Depiction of the basics of EEG and MEG. (a) An ionic current generated bylargeneuronalcells. Theprimarycurrentisrepresentedbyablacklargearrowwhile the dark lines represent secondary currents. Magnetic fields which are generated by primary and secondary currents are represented by dashed lines. (b) The neuronal pop- ulations generate mass effect of currents which behave as a current dipole (red arrow). EEG detects the secondary currents (circular yellow lines) induced by primary genera- tor. MEG captures the magnetic fields which travel more easily within tissues. From [Bai12]. of SQUIDs (Superconducting Quantum Interference Devices), MEG is able to detect magneticfieldsontheorderofseveraltensoffemtoTeslas[Koc08]. AtypicalMEGset upisshowninFig. 1.12. (a) (b) Figure 1.11: (a) An elastic EEG cap with 60 electrodes. From [Bai12]. (b) EEG elec- toresrelativetothehumanheadandthefoldedcortex. From[Nun95]. 18 Figure1.12: AtypicalMEGequipment. MEGmeasurementsareobtainedinsideamag- netically shielded room (left). Liquid helium (upper right) stored in a Dewar (left and upper right) keeps MEG sensors which use low-electronic temperature cool. Visualiza- tionoftime-evolvingmagneticfieldtopographies(lowerright). From[BML01]. In addition to the temporal information, spatial information is also required to under- standbrainactivity. Unfortunately,EEGandMEGofferverylimitedspatialresolution. One of the reasons of this limitation is the relatively small number of spatial measure- ments (64 256 electrodes for EEG and 306 sensors for MEG). Spatial resolution can be enhanced by mapping the EEG/MEG measurements on the cortex. Advanced sig- nal processing approaches have been proposed to reconstruct the sources [BML01]. This mapping involves two separate but related problems: the forward problem and theinverseproblem. 19 The forward problem determines the scalp potentials and external fields using the quasi-static approximations of Maxwell’s equations for a given source distribution and head volume conduction properties. It is a well-defined problem with a unique solu- tion for both EEG and MEG [Plo69, MP95]. Since this is a linear mapping, it can be representedbyatransfermatrixalsocalledleadfieldmatrix. Thegoalofinversemodelingistoreconstructthedistributionofthecurrentsources within the brain given the EEG/MEG measurements and head volume conductor prop- erties. This can also be viewed as mapping from sensor space to source space. Due to thelimitationofthenumberofelectrodesorsensorsinEEG/MEG,thisproblemisfun- damentallyill-posed[FET00,GDMH + 00]. Inotherwords,thereareaninfinitenumber of solutions that fit the data well. Uniqueness can be achieved by incorporating a priori information such as anatomic constraints on the sources [DS93], on volume conductor [HMO + 87,HS89],orfunctionalconstraintsprovidedbyotherimagingmodalitiessuch asfMRI [DLF + 00, LH08, HL08]. Although inversemodeling improvesthespatial res- olutionofEEG/MEG,theambiguityofthesolutionisstillpresent. Thismakesstudying EEG/MEGrecordingschallenging. Inthisdissertation,wewilldescribeapproachesthat areproposedtobeappliedtoEEGandMEGrecordings. 1.2.5 ComparisonofImagingTechniques The functional neuroimaging techniques described here can be classified into two cate- gories: (1)metabolicandhemodynamic(2)electro-magnetic[HTM99]. PETandfMRI belongs to the former, whereas ECoG, LFP, single units, EEG and MEG belong the electro-magnetic techniques. These two categories measure fundamentally different mechanisms of the neurons despite the relationship between them. Furthermore, each techniquehasdifferentdrawbacksrelatingtospatialandtemporalresolution[CFM + 10]. fMRI has dominated neuroimaging studies owing to the excellent spatial resolution it 20 offers [HSM04]. However, it is an indirect measure of neruonal activity and the tem- poral resolution is poor for capturing the neuronal changes which occur in the order of milliseconds. ECoG, LFP and single unit measurements are rich in spatial and tempo- ral resolutions but they require invasive surgery [GG01], hence they are not applicable to the healthy subjects. MEG and EEG are a direct measure of neuronal activity and are able to capture temporal changes in the neuronal activity despite the limited spatial resolution[BML01]. In this dissertation we focus on the elctro-magnetic imaging techniques to develop newsignalprocessingandstatisticaltoolsforfunctionalconnectivity. 1.3 ConnectivityintheBrain Structuralsegregationandthefunctionalspecializationandintegrationinthebrainhave been well-established for centuries [Fri11]. Different neuronal populations distributed throughout the cerebral cortex in the brain are specialized in processing different neu- ronal computations [TE98, Hor03]. Processing a single function may involve coordi- nation between different specialized regions [Fri09]. Brain connectivity networks are thought to reveal the physiological basis for segregation and integration in the brain [GBR08]. Detectingandquantifyingtheinteractionsbetweendistinctneuronalpopula- tions can lead to important insights into the large scale complex networks that underlie highcognitiveprocesses. Thisisalsocrucialinunderstandingtheclinicalconsequences ofsomeofthebraindisorders[GBR08]. Two types of brain connectivity have been defined: structural connectivity and func- tionalconnectivity. Mappingthebrainnetworksusingdifferentmodalitiesisalsocalled 21 Figure 1.13: Structural and Functional Brain Networks. First, the network nodes are defined for both connectivity types. Second, measures of interaction both for structural and functional connectivity are computed. This could be the correlation metrics based on temporal interactions such as coherence, Granger causality for functional networks while the connection probability between two regions of an individual diffusion tensor imaging data set, or the inter-regional correlations in cortical thickness or volume MRI measurements estimated in groups of subjects for structural networks. Finally, connec- tivitynetworksareobtainedafterstatisticalthresholding. From[BS09]. human connectomics [STK05]. Fig. 1.13 illustrates the steps to estimate both func- tional and structural connectivity networks. We describe the structural and functional connectivitynetworksindetailbelow. 22 1.3.1 StructuralConnectivity Structural networks analysis have provided insight into the physical or structural path- waysorlinksbetweenbrainregions. Thecharacterizationofthelinksaredeterminedby statisticalmeasuresofassociationsincludinginterconnectingaxonalpathways;orinter- regionalcross-correlationinanatomicalstructuressuchascorticalgraymatterthickness graymattervolume,andsurfacearea[BS09]. Theformercanbeobtainedbyexamining the white matter fiber connections among gray matter regions from diffusion MRI data [HE10]. Thelatterone,ontheotherhand,usesstructuralMRIdata. StructuralMRI The rich morphological information of structural MRI has been used to charaterize anatomical connectivity of the brain. Various studies have reported correlation in gray mattermorphologybetweenanatomicallyandfunctionallylinkedbrainregions[HE10]. He et al. [HCE07] used the statistical correlation of cortical thickness measure- mentsin124individualstoinvestigatethenetworkofanatomicalconnections. Cortical thickness was chosen as a morphometric measure as it comprises the structural infor- mation of the cells [TC69, NBT + 05]. Significant correlation values are considered to be an indicator of the presence of corticocortical pathways. Recently, Sanabria-Diaz et al. [SDMGIM + 10] used surface area in addition to the cortical thickness to build net- works of structural connectivity and demonstrated distinct properties of the interaction betweenbrainregions. DiffusionMRI Structural networks in the brain have been mapped using recent noninvasive imaging methodsdiffusionMRIandtractography[HE10]. Diffusiontensorimaging(DTI)quan- tifies diffusion of free water molecules in three dimensions in biological tissue. The 23 motion of molecules are supressed by the architectural properties of tissues such as axonsinthebrain. Itisknownthatwaterdiffusionishighlyanisotropicincertainwhite matterregions,withpreferredmovementalongnervefibers[IMCRMG + 07]. HenceDTI reflects the arrangement and orientation of the fiber pathways [GBR08]. The anatomi- cal circutry of the living brain that DTI provides is basically the white matter pathways betweenvoxels. Several studies have demonstrated that DTI can produce structural network results consistent with established pathways formed by major tracts in the brain [CLC + 99, LB03, CJD + 03] despite some limitations in data acquisition and processing algorithms [DBSB + 05]. 1.3.2 FunctionalConnectivity Quantifying dynamic brain networks provide greater insights into the understanding of the brain. Modalities introduced in section 1.2 provide characterization of activity in the brain in the form of time series. Functional connectivity analysis seek to elucidate temporal interaction between these time series from functionally segregated regions. A multivariatenetworkcanbeconstructedbycomputingtheinteractionmeasuresbetween assignednodesinthebrain. Functional connectivity studies can be categorized in two distinct types. One type aimstofindfunctionalconnectivitynetworksunderanexplicittask. Thesestudiesoften require aggregated data over 10 or more subjects which might not be possible in some cases such as clinical applications. Furthermore, most of the patients may not be able to perform some particular tasks efficiently. Therefore, it is difficult to make general- izationandinterpretationfromthistypeofconnectivitynetworks. Despitetheselimita- tions, functional connectivity studies have focused on the brain’s response to a task or stimulihavedominatedneuroscienceforthepast50years. However,thebraindoesnot 24 stop working even in the absence of external input [FR07]. The motivation behind our interest in resting state can be explained through the energy consumption of the brain. It is known that the resting human brain consumes 20% of the body’s energy and this value increases only by 5 % in response to task which is relatively small compared to the resting energy consumption [RM06]. This type of connectivity is so-called “rest- ingstate”connectivityandhasreceivedsignificantattentioninrecentyears[HWA + 13]. Understandingthecomponentofthebrainwhichconsumesmostofeneryiscrucialfor acomprehensiveunderstandingofbraindynamics. Bloodoxygenationleveldependent(BOLD)fMRIhasattractedmanyresearchersas itprovidesnon-invasiverepresentationofthebrainactivitywithhighspatialresolution. Biswaletal. [BZYHH95]demonstratedcorrelationbetweenlowfrequencyfluctuations of functionally related brain regions in the resting state. This study introduced a new approach to neuroscience by showing that a resting state activity is not random. Their results have been replicated by several groups [CHA + 00, LMS98, CHA + 01, XPGF99, FSZR05, DLSDS + 05]. Fox et al. [FCS + 06] discovered the connectivity pattern using correlationbetweensimultaneoussignalsintheabsenceoftaskwhichwasalsoobserved inresponsetoexternaltasksorstimuli. As discussed in section 1.2.5, BOLD fMRI measurements are indirect assessment ofbrainactivity. Theyarelimitedbypoortemporalresolutionsincetheyrelatetoblood flow. The frequency range of these measurements is upper bounded by 0.1 Hz. How- ever, electrophysiological signals in the brain have a wide range of frequency spectrum which are traditionally divided into different bands (delta, 1-4 Hz; theta, 5-8 Hz; alpha, 9-14Hz;beta,15-24Hz;gamma,>24Hz). Electromagneticmeasurementsofthebrain signals using single unit, LFP, ECoG, EEG/MEG allow us to study brain activity in a wider frequency range due their ability to capture temporal variations in the range of miliseconds. Electrophysiologicalsignalscanbeusefullycharacterizedintermsoftheir 25 oscillatory components either through band-pass filtering into the standard frequency bands or using broadband spectral representations of the data. Interactions can then be analyzedusingmeasuresofwithinandbetweenfrequency-bandcouplingbetweenelec- trodeormagnetometerpairs. InvasivemodalitiessuchasLFPandECoGhavealsoarich spatialresolution. Furthermore,connectivityarchitecturecharacterizedbythesemodal- itiesareshowntobelinkedtostimulus-evokedBOLDresponsesandrestingstatefMRI [LPA + 01, HSZ + 08, SL08, NMD + 08]. Although invasive recordings provide such high temporal and spatial resolution, they require surgery and therefore not ideal for healthy subjects. Not only do these methods require invasive surgery, but the measurement devicescancoveronlyasmallfractionofthecorticalsurface. Noninvasive techniques such as EEG/MEG allow us to measure electromagnetic activity over most of the cortex through the electrodes placed outside the head. Rest- ing state networks can be characterized using these modalities in sensor domain [BMLA + 06,HRPH10,LLG06,LFdZD10,MPDG + 07,MBM + 10]. Despitethewhole- head coverage and excellent temporal resolution, these modalities have poor spatial resolutions. If EEG or MEG data are first mapped back onto cortex using an inverse mapping procedure [BML01], then we can improve the spatial resolution and com- pute interactions between time series averaged over cortical regions of interests (ROIs) [HBB + 12]. In order to quantifty functional connectivity, several approaches have been pro- posed including linear measures such as cross-correlation, coherence [Car87]; nonlin- ear measures such as mutual information [Rou99] and phase locking value [LRM + 99, TRW + 98] and directional measures such as Granger Causality [BDL + 04], Directed Transfer Function [KDTB01] ,and Partial Directed Coherence [KKB04]. This abun- danceoftheavailabletechniquesforfunctionalconnectivityispromisingforthefuture ofneuroimagingbutalsoconfusing. 26 In this dissertation we restrict our attention to coherence and phase locking value computed between pairs of electrodes, magnetometers or cortical ROIs. Coherence is a complex measure of phase and amplitude similarity computed as a function of fre- quency [NSW + 97, KSJS06, CK91]. Phase locking value considers only the relative phasebetweenthetwosignals[TRW + 98]. Inthisdissertation,wewillclarifythestatis- ticalrelationshipsbetweentheseinteractionmeasures. Wewilllateraddresstheproblem ofcoherenceinpresenceofcross-talkbetweensignals. Due to the limited spatial resolution of EEG and MEG, it is likely that the coher- ence between electrodes or sensors do not reflect the underlying network. Neighbor- ing electrodes or sensors might pick up the common source within the brain resulting spuriosinteractionbetweenthesetimeseries. Thisproblemistermed“volumeconduc- tion” or “cross-talk” problem [NSS + 99] and have been increasingly addressed recently [GHAEC08, GMDV + 11, HPBL10, MDGN08, NZN + 08, NZN + 08, SG09, SND07]. Similarly,commonreferenceelectrodeinEEGcontributesasignificantamountofinter- action[FRBM88]. Thus, the resultswillbeaffectedbyrelativepowerchanges. Inverse modelingofEEGandMEGdoesnotremovethecross-talkentirely[SP98]. Ourgoalistoreviewanddevelopasolidmathematicalunderstandingtotheexisting measures which are proposed to overcome the problem of volume conduction. We will proposeanovelapproachwhichismorerobusttocross-talk. 1.4 BrainDisorders Thissectionreviewssomeofthebraindisordersthatshowedalteredfunctionalconnec- tivitynetworksbasedonthepreviousstudies. 27 1.4.1 Schizophrenia Significantnumberofstudieshavereportedirregularfunctionalconnectivityofpatients with schizophrenia [HWA + 13, Gre08]. Recently, Sakoglu et al. performed func- tional connectivity analysis of spatially independent networks in healthy controls and schizophrenia patients during an auditory oddball task [SPK + 10]. They showed group differences in time-frequency patterns of networks between two groups. Furthermore, they reported that sensory, motor, and frontal networks had less engagement with the othernetworksconsistentwiththeexistinghypothesisthatschizophreniapatientsshow lesssegregatedmotor,sensory,cognitivefunctionsandlesssegregateddefaultmodenet- work activity when engaged with a task. Bassett et al. reported a relationship between univariate and bivariate measures which has been significantly altered in the patients withschizophrenia[BNM + 12]. 1.4.2 MajorDepression Severalstudieshavesupportedthealterationinfunctionalconnectivityinpatientsdiag- nosedwithmajordepressionbutmoreworkisneededtobetterunderstandthedynamics of this disorder [HWA + 13]. Horn et al. reported the relationship between increased functionalconnectivitybetweentwobrainregions(pgACCandleftanteriorinsularcor- tex (AI)) and glutamatergic levels in pgACC [HYS + 10]. This relationship was more significant in severely depressed patients while was not observed in healthy controls. Recently,Hamiltonetal. utilizedGrangercausalitytoidentifyasetofstructureswhich indicatesdifferentcausalitynetworkwithvACCindepressedpatientsthanhealthysub- jects [HCT + 10]. Allen and Cohen reported the reduction of alpha synchrony in resting stateEEG[AC10]. 28 1.4.3 Alzheimer’sdisease Various studies have reported altered functional connectivity networks in Alzheimer’s patients compared to healthy controls [HWA + 13]. Recently, Jones et al. identified the temporal differences of functional conectivity between patients with Alzheimer’s dis- ease and age-matched healthy individuals which is explained by non-stationary nature ofthebrainsmodularorganization[JVM + 12]. Thiswasone ofthefirststudiesdemon- stratingRS-fMRIchangesinAlzheimer’spatientsbeyondaverageFCmetrics. Previous studies reported altered static functional connectivity networks in Alzheimer’s patients comparedtohealthysubjects[BST + 09,SMR + 08,SRM + 07]. 1.4.4 Epilepsy Epilepsy is a notable example for brain disorders with a devastating feature of seizure. The seizure have been characterized by the interaction between signals recorded from distinct sites at seizure onset [BWV + 99, BWR + 04, BWB + 01, DS92, GL96]. Various studies have elucidated functional connectivity disturbances indicating epileptic abnor- malities[ETS + 13]. Therefore,functionalconnectivityisanimportantphenomenonthat can be used to identify epileptic networks. Reliable estimation of these networks will help treatment planning for patients. Engel et al. provided a comprehensive review on connectomics at the whole brain including functional studies including the various approaches to resting fMRI, invasive, and noninvasive EEG as well as structural con- nectivitystudies[ETS + 13]. 29 1.5 OrganizationoftheDissertation The organization of the dissertation is as follows: we describe the first part of our work in chapter 2: phase synchronization. We discuss the properties of the phase synchro- nization measures including phase locking value, phase lag index and coherence. After showing one to one relationship between phase locking value and coherence for Gaus- siansignals,wewillthenturntothemultivariateextensionofPLVandshowitsrelation to the partial coherence in chapter 3. In chapter 4 we discuss the cross-talk problem in electrophysiological signals. We review the existing measures including lagged coher- ence, phase lag index, imaginary coherence and then propose a novel measure “partial lagged coherence (PLC)” as a robust measure against multiple cross-talk. We conclude by proposing future work including the development of partial phase locking value and applicationtoelectrophysiologicalrecordings. 30 Chapter2 PhaseSynchronization 2.1 Introduction Phase locking is a fundamental concept in dynamical systems that has been used in control systems (the phase-locked loop) and in the analysis of nonlinear, chaotic and nonstationary systems. Since the brain is a nonlinear dynamical system, phase locking is an appropriate approach to quantifying interaction. A more pragmatic argument for its use in studies of LFPs (Local Field Potentials), EEG and MEG is that it is robust to fluctuationsinamplitudethatmaycontainlessinformationaboutinteractionsthandoes therelativephase[LRM + 99,MLDEE00]. The most commonly used phase interaction measure is the Phase Locking Value (PLV), the absolute value of the mean phase difference between the two signals expressedasacomplexunit-lengthvector[LRM + 99,MLDEE00]. Ifthemarginaldistri- butionsforthetwosignalsareuniformandthesignalsareindependentthentherelative phase will also havea uniform distributionand the PLVwill be zero. Conversely, if the phases of the two signals are strongly coupled then the PLV will approach unity. For event-related studies we would expect the marginal to be uniform across trials unless thephaseislockedtoastimulus. Inthatcase,wemayhavenonuniformmarginalwhich couldinprincipleleadtofalseindicationsofphaselocking. Whencomparingelectrodepairsthatshareacommonreferenceoroverlappinglead field sensitivities, or when investigating cortical current density maps of limited reso- lution, the PLV suffers from sensitivity to linear mixing in which the same source can 31 contribute to both channels. In these cases, the PLV can indicate an apparent phase locking with the relative phases concentrated around zero. [SND07] proposed an alter- nativemeasure,thePhaseLagIndex(PLI)thatisrobusttothecommonsourceproblem. PLI quantifies the asymmetry of the relative phase distribution about zero and so will producelargevaluesonlywhentherelativephaseispeakedawayfromzero. In this dissertation we first define nonparametric estimates of PLV and PLI, and consider the bias intrinsic in the sample PLV estimator. We derive an expression for the unbiased estimator of the squared PLV and show equivalence to the Pairwise Phase Consistency (PPC) metric recently proposed by [VVWW + 10]. We then investigate the relationship between PLV and PLI and two possible parametric distributions that can be used to model relative phase. The first of these, the von Mises distributions, is the maximum entropy distribution over the class of circular distributions [JS01]. The sec- ondmodelistherelativephasedistributionassociatedwithcomplexcircularlysymmet- ric Gaussian processes. This model is appropriate for complex signals generated fom jointly Gaussian real signals through use of the Hilbert transform. The relative phase distribution is obtained by marginalizing the joint Gaussian distribution with respect to the amplitude of the two complex signals. We derive closed-form expressions for the relationshipbetweenPLVandtheparametersofthevonMisesandGaussianmodels. Invasivemicroelectroderecordingscanbeusedtoinvestigatebothmultiunitactivity, which reflects axonal firing rates, and the local field potentials (LFPs) associated with dendritic and volume conduction currents. In this paper we are concerned with the application of PLV and PLI measures to LFPs as well as noninvasive EEG and MEG measurements that similarly result from dendritic and volume conduction currents. We use LFP recordings from a macaque monkey study [BDY99] to investigate whether the von Mises and Gaussian distributions are appropriate for modeling relative phase between pairs of electrodes. We then compare the ability of two different estimators of 32 PLV, associated respectively with the von Mises and Gaussian models, to detect phase lockingbetweenelectrodes. The goal of this chapter is to clarify the relationships between nonparametric esti- matorsofPLVandPLIandtwowell-knownparametricdistributionsthatcouldbeused to model phase interactions. A second goal is to investigate the relationship between PLV and coherence (cross-correlation for complex narrow-band signals) when analyz- ing LFP data. We beginby stating, and where appropriate deriving, these relationships. WethenpresentcomputationalsimulationsandanalysisofexperimentalLFPdatausing differentPLVestimators. 2.2 MeasuresofPhaseSynchronization 2.2.1 ThePhaseLockingValueandPhaseLagIndex Phase synchronization between two narrow-band signals is frequently characterized by the Phase Locking Value (PLV). Consider a pair of real signals s 1 (t) and s 2 (t), that have been band-pass filtered to a frequency range of interest. Analytic signals z i (t) = A i (t)e jϕ i (t) for i = f1;2g and j = p 1 are obtained from s i (t) using the Hilbert transform: z i (t) =s i (t)+jHT(s i (t)) (2.1) whereHT(s i (t))istheHilberttransformofs i (t)definedas HT(s i (t)) = 1 P.V. ∫ 1 1 s i (t) t d (2.2) and P.V. denotes Cauchy principal value. Once the analytic signals are defined, the relativephasecanbecomputedas 33 ∆ϕ(t) = arg ( z 1 (t)z 2 (t) jz 1 (t)jjz 2 (t)j ) : (2.3) TheinstantaneousPLVisthendefinedas[LRM + 99,Cel07] PLV(t)≜ E [ e j∆ϕ(t) ] (2.4) where E[:] denotes the expected value. The PLV takes values on [0;1] with 0 reflecting the case where there is no phase synchrony and 1 where the relative phase between the two signals is identical in all trials. PLV can therefore be viewed as a measure of trial to trial variability in the relative phases of two signals. In this work we use the Hilbert transform but the continuous Morlet wavelet transform can also be used to compute complex signals, producing separate band-pass signals for each scaling of the wavelet. [QKKG02]and[LVQFL + 01]haveshownthatbothapproachesyieldsimilarresults. When computing synchrony between pairs of electrodes or cortical locations, nonzero PLVs can arise from a single source contributing to both signals as a result of either volume conduction in channel space or limited spatial resolution in the case of cortical current density maps [NSW + 97, TRW + 98, DGCV02, GVN + 05, ARN + 05, VOvW + 11]. In this case of direct linear mixing there is no phase lag between the two signals potentially resulting in a large value of PLV. Linear mixing can therefore easily be mistaken for phase locking between distinct signals. To distinguish these two con- ditions we need a different measure of phase locking that is zero in the case of linear mixing but nonzero when there is a consistent nonzero phase difference between the two signals. The Phase Lag Index (PLI) [SND07] achieves this goal by quantifying the asymmetryofthedistributionofrelativephasearoundzeroandisdefinedas PLI≜jE[sign(∆ϕ(t))]j: (2.5) 34 PLI takes values on the interval [0;1] and is zero if the distribution of relative phase is symmetricabout 0or. InpracticePLVandPLIaretypicallyestimatedbyaveragingovertrialsand/ortime [LRM + 99, LRM + 00, MLDEE00, SND07, ABES10]. For notational convenience, we will drop the explicit dependence on t in the following. A nonparametric estimate of PLVcanbecomputedbyapproximatingequ. (2.4)byaveragingovertrials: ^ PLV sample ≜ 1 N N ∑ n=1 e j∆ϕn(t) (2.6) where n indexes the trial number and N is the total number of trials. The estimator generalizesinanobviouswaytoincorporateaveragingovermultipletimesamples. The correspondingnonparametricestimatorforPLIis ^ PLI sample ≜ 1 N N ∑ n=1 sign(∆ϕ n (t)) : (2.7) In the following section we consider the relationship between PLV and PLI and the parameters of two alternative probability distributions that can be used to characterize phase interactions: the von Mises and the bivariate circularly symmetric Gaussian. We firstconsidertheissueofbiasinthenonparametricPLVestimator. Without specifying the distribution of relative phase, we cannot find an expression forthebiasinthesamplePLVdefinedin(2.4). However,asshownby[VVWW + 10],in thegeneralcase ^ PLV sample isabiasedbutconsistentestimatorofPLV: E [ ^ PLV sample ] > PLV (2.8) lim N!1 ^ PLV sample = PLV: (2.9) 35 Working with the squared estimator of PLV, rather than PLV itself, it is straightforward to derive the following nonparametric expression for the bias of ^ PLV 2 sample as a func- tionofN (seeAppendixA.1): E [ ^ PLV 2 sample ] = 1 N + ( 1 1 N ) PLV 2 : (2.10) This expression holds regardless of the distribution of the relative phase. Rearranging theexpressionweobtainthefollowingunbiasedPLV 2 estimator: ^ PLV 2 ub sample ≜ 1 N1 ( ^ PLV 2 sample N1 ) : (2.11) Interestingly, this measure is identical to Vinck’s Pairwise Phase Consistency (PPC) measure PPC = 2 N(N1) N1 ∑ n=1 N ∑ m=n+1 cos(∆ϕ m (t)∆ϕ n (t)) (2.12) as we show in Appendix A.2. While this is an unbiased estimator, in general its square root is not an unbiased estimator of PLV. The bias in ^ PLV sample is dependent on dis- tribution and we were unable to find closed form unbiased estimators for either the von MisesorcircularGaussiandistributionsthatweinvestigatebelow. We illustrate the bias and variance of these estimators in Fig. 2.1 for relative phase values sampled from the von Mises distribution (which we describe in detail in the followingsection). ByvaryingtheconcentrationparameterofthevonMisesdistribution we are able to produce differing values of PLV. The figure shows that bias decreases rapidly with number of samples and is negligible for N > 50. We also see that bias reduces as true PLV increases for fixed N. Finally, as observed by [VVWW + 10], we seethatthereisasmallincreaseinvariancewhenusingtheunbiasedratherthanbiased 36 (a) (b) Figure 2.1: Bias and variance of ^ PLV 2 sample (red) and ^ PLV 2 ub sample (blue) as a func- tion ofN, the number of samples used to compute the estimators: (a) mean value vsN forfourdifferenttruevaluesofPLV;(b)variancevs. N forthesametruevaluesofPLV. Samples were drawn independently from the von Mises distribution for four different concentrationparametervalues. measure. Sincebiasissmallexceptwhenthenumberoftrialsissmall,inthefollowing wewillcontinuetoworkwiththebiasedestimatorofPLV,whichsimplifiesouranalysis of the relationships between PLV and the parameters of the von Mises and circular Gaussianmodels. 2.2.2 PhaseLockingValueandthevonMisesDistribution The von Mises distribution is the most widely used model for circular (or periodic) random variables, in part because it is the circular distribution with maximum entropy subjecttoconstraintsonitsfirsttrigonometricmoments[JS01]. Theprobabilitydensity function(pdf)ofthevonMisesdisributionis: p(ϕj;) = 1 2I 0 () e cos(ϕ) (2.13) 37 with concentration parameter2 [0;1) and mean2 [ ] ;I 0 () is the modi- fiedBesselfunctionofzerothorder I 0 () = 1 2 ∫ 2 0 e cos() d: (2.14) The von Mises distribution is unimodal and symmetric about as illustrated in Fig. 2.2(a) for = 0. As the concentration parameter increases, the relative phase ϕ becomesincreasinglyconcentratedaboutitsmean. Inthelimitasgoesto0,itreduces to a uniform distribution corresponding to the case where the phases of the two signals aremutuallyindependent. (a) (b) Figure 2.2: Probability density functions: (a) von Mises distribution with mean = 0 for a range of concentration values, ; (b) relative phase distribution for bivariate cir- cularly symmetric Gaussian models for different values of cross-correlation magnitude jR 12 j = 12 p 11 22 withphase 12 = 0. ThePLVforthevonMisespdfcanbefoundfromthemomentgeneratingfunction: PLV vonmises = E [ e j∆ϕ(t) ] = I 1 () I 0 () (2.15) 38 as shown by [JS01]. The PLV is therefore a monotonic function of and independent of. Incontrast,PLIcanbefoundas: PLI = E[sign(∆ϕ(t))] = ∫ 0 p(ϕj;)dϕ+ ∫ 0 p(ϕj;)dϕ : (2.16) (a) (b) Figure2.3: Plotsofsample(a)PLVand(b)PLIasafunctionofconcentrationparameter and mean for samples drawn from the von Mises distribution. Note that while PLV isindependentof,PLIdependsonbothand. A closed form for this expression is intractable but it is clear that unlike PLV, PLI is a functionofbothand. WeshowthisdependenceinFig. 2.3forsamplesdrawnfrom the von Mises distribution. This figure illustrates the monotonic relationship between PLVandandthemorecomplexinteractionbetween (;)andPLI. From equ. (2.15) we see that the PLV can be computed from an estimate of the concentration parameter . As shown in [JS01], the maximum likelihood estimator of PLVforthevonMisesdistributioncanbecomputeddirectlyfromamaximumlikelihood 39 estimator of and furthermore, the resulting estimator is identical to the sample PLV estimatorinequ. (2.6): ^ PLV vonmises ≜ I 1 ( ML ) I 0 ( ML ) = 1 N N ∑ n=1 e j∆ϕn(t) = ^ PLV sample : (2.17) 2.2.3 Phase Locking Value and Circularly Symmetric Gaussian Processes We now turn to an alternative statistical model for relative phase. Under the assump- tionthatthetimeserieswhosephasesarebeingcomparedarejointlyGaussian,thenthe distribution of their relative phase, and hence the PLV and PLI, must be governed by thepropertiesandparametersofthatGaussianprocess. Wenowexaminetheserelation- ships. Let the sources s 1 (t) and s 2 (t) be jointly Gaussian zero mean processes. Then, the complex random vectorz(t) = [ z 1 (t) z 2 (t) ] T wherez i (t) = s i (t) +jHT(s i (t)) = A i (t)e jϕ i (t) follows a circularly symmetric complex Gaussian distribution that satisfies the condition E [ z(t)z(t) T ] = 0 and has the property that the real and imaginary com- ponents ofz(t) are mutually independent [Gal08]. The pdf of the circularly symmetric complexGaussiandistributionforavectorz(t)is p(z(t)) = 1 2 jK z j exp { z(t) H K 1 z z(t) } (2.18) where, for the bivariate caseK z △ = 2 4 11 12 e j 12 12 e j 12 22 3 5 1 = E [ z(t)z(t) H ] is the covariancematrix ofz(t), with 0 12 <1 and 12 . Inverting this matrix and normalizing gives the coherence (cross-correlation in this case) between z 1 (t) and z 2 (t)asR 12 ≜ 12 p 11 22 e j 12 . 40 We rewrite the pdf in polar coordinates in order to investigate the relative phase distribution: p(A(t);(t)) = A 1 (t)A 2 (t) 2 jKzj exp { [ A 1 (t)e jϕ 1 (t) A 2 (t)e jϕ 2 (t) ] 2 4 11 12 e j 12 12 e j 12 22 3 5 2 4 A 1 (t)e jϕ 1 (t) A 2 (t)e jϕ 2 (t) 3 5 9 = ; = A 1 (t)A 2 (t) 2 jKzj exp { 2 ∑ i=1 ii A i (t) 2 2 12 A 1 (t)A 2 (t)cos(ϕ 1 (t)ϕ 2 (t) 12 ) } (2.19) whereA(t)≜ [ A 1 (t) A 2 (t) ] and(t)≜ [ ϕ 1 (t) ϕ 2 (t) ] . For the univariate case the circularly symmetric Gaussian pdf in polar form is separable. In other words, if we write z = Ae iϕ , then the density can be writ- ten p(z) = p(A;ϕ) = p(A)p(ϕ), where p(A) follows a Rayleigh distribution and p(ϕ) is uniform [DR58]. However, this is not the case for the bivariate case, i.e. p(A(t);(t)) ̸= p(A(t))p((t)) , but rather the phase and amplitude are mutually coupled unless the cross-correlation is 0. We therefore now consider two different dis- tributions for the relative phase ∆ϕ(t) = (ϕ 1 (t)ϕ 2 (t)): the distribution conditioned ontheamplitudes,andthedistributioninwhichwemarginalizeouttheamplitudes. FromBayesrule,thephasedistributionconditionedonamplitudecanbewrittenas: p((t)jA(t)) = p(A(t);(t)) p(A(t)) : (2.20) As shown in Appendix A.3, this results in the following conditional distribution for relativephase: p(∆ϕ(t)jA(t)) = 1 2I 0 (2 12 A 1 (t)A 2 (t)) e 2 12 A 1 (t)A 2 (t)cos(∆ϕ(t) 12 ) : (2.21) Note that this has the form of a von Mises distribution, but with the mean and concen- tration parameters a function of not only the parameters of the Gaussian but also the amplitudeoftheobservations. Thisresultclearlyshowsthatwhilewemaybeinterested 41 only in the relative phases of two signals, amplitude and phase are not independent, so that the amplitudes of the signals affect the relative phase distribution (and vice versa). Put another way, the amplitude and phase of jointly Gaussian signals do not contain independent information about those signals. Because of this amplitude dependence, to determine the distribution of the relative phase alone we need to marginalize with respect to amplitude, as also independently shown by [VON11], to obtain the relative phasedistribution: p(∆ϕ(t)) = 1jR 12 j 2 2(1 2 ) ( 1 cos 1 ( ) √ (1 2 ) ) (2.22) whereR 12 ≜ 12 p 11 22 e j 12 and ≜jR 12 jcos(∆ϕ(t) 12 ). Notethatinthiscasethe relative phase is solely a function of the cross-correlationR 12 betweenz 1 (t) andz 2 (t), whichalsogovernsthecorrelationinamplitudebetweenthetwosignals. Fig. 2.2(b)showsplotsoftherelativephasedistributionforthecircularlysymmetric complex Gaussian model for different values of R 12 with 12 = 0. Varying 12 from zerowillcircularlyshiftthedistributionsothatitsmaximumisat 12 . Incomparisonto the von Mises distribution,jR 12 j plays the role of the concentration parameter, and 12 the mean. However, there are differences in the shapes of these two distributions, with theGaussiancaseshowingslightlylongertailedbehavior. Usingthemomentgenerating function of the relative phase PDF (equ. 2.22), as we show in Appendix A.4, we find thephaselockingvalueforthecircularlysymmetricGaussianmodelis PLV circgauss △ = E [ e j∆ϕ(t) ] = p 2 ( 1 11 22 2 12 )[ w 3=2 2 F 1 ( 3 4 ; 5 4 ;1;w 2 ) + 3 4 w 5=2 2 F 1 ( 5 4 ; 7 4 ;2;w 2 )] (2.23) 42 Figure2.4: PlotofthemonotonicrelationshipbetweenPLVandthemagnitudeofcross- correlation,jR 12 j = 12 p 11 22 ,forthebivariatecircularlysymmetricGaussianmodel. where w ≜ 2 12 2 11 22 2 12 and 2 F 1 (:) represents the hypergeometric function. Note that PLV is a function of the magnitude of the cross-correlationjR 12 j = 12 = p 11 22 . Fig. 2.4 shows this one-to-one relationship. The significance of this result is as follows. The parameter 12 = 11 22 represents the magnitude of the cross-correlation between the two complex Gaussian processes z 1 (t) and z 2 (t). If these are narrow-band signals, the cross-correlation at a single frequency is equivalent to coherence between the two processesatthatfrequency. Consequently,PLVandcoherenceareequivalentmeasures. Of course, this holds only for the Gaussian case, and we will return to the question of whether electrophysiological signals can be adequately modeled as jointly Gaussian in thefollowingsection. As with the von Mises case, we do not have a closed form for PLI using the integral in equ. (2.16) after substituting in the relative phase distribution in equ. (2.22) for the Gaussian case. However, we can determine the relationship between PLI and the Gaussianparameters(jR 12 j; 12 )byMonteCarlosampling. Wegeneratedsamplesfrom 43 (a) (b) Figure2.5: Plotsofsample(a)PLVand(b)PLIasafunctionofcross-correlationmag- nitudejR 12 j = 12 p 11 22 and phase 12 for the bivariate circularly symmetric Gaussian distribution. SimilarlytothevonMisesdistribution,PLVisindependentofphasewhile PLIdependsonbothcross-correlationmagnitudeandphase. Samplesweredrawnfrom relativephasedistributioninequ. (2.22). the relative phase distribution in equ. (2.22) for a range of values ofjR 12 j and 12 . Fig. 2.5 shows sample estimates of PLV and PLI as a function of these two parameters. Again,similarlytothevonMisescase,whilePLVdependsonlyonthejR 12 jvalue,PLI isafunctionofbothjR 12 jand 12 . The relationship in equ. (3.6) between coherence and PLV also gives us an alterna- tiveestimatorforthelatter. RatherthandirectlycomputingsamplePLV,wecaninstead compute the sample estimator ofjR 12 j and then substitute this into equ. (3.6). We refer to this estimator as ^ PLV circgauss . In the following section we compare this with the samplePLVestimateinbothsimulatedandexperimentalLFPdata. 44 2.3 Results 2.3.1 SimulationsBasedontheGaussiandistribution IntheprevioussectionwedescribedthreedifferentestimatorsforthePLV: ^ PLV sample (equ. 2.6), √ ^ PLV 2 ub sample (equ. 2.11), and ^ PLV circgauss (the PLV computed from the sample coherence, Equ. 3.6). Their relative performance in terms of bias and vari- ance will depend on the true distribution of the data. Before considering the case of experimental electrophysiological data, we will first explore their behavior for simu- latedbivariateGaussianprocesses. (a) (b) Figure 2.6: Plot of bias and variance for the sample PLV, root of unbiased sample PLV, and PLV computed from sample coherence. (a) plot of the mean over 1,000 Monte Carlo trials as a function of number of samples used to estimate the parameter for two different values of coherence; (b) corresponding plot of variance for each of the three estimators for the two different coherence values. Samples were drawn independently frombivariateGaussianprocesses. WegeneratedindependentsamplesfromtwodifferentbivariateGaussianprocesseswith normalized coherence values of jR 12 j = 0:25 and jR 12 j = 0:91, with analytic PLV 45 values of 0:20 and 0:83, respectively. This was achieved by multiplying pairs of inde- pendent zero mean uncorrelated Gaussian random variables by the inverse square root of the 2 x 2 covariance matrix 2 4 1 R 12 R 12 1 3 5 . We then generated the complex signal by applying the Hilbert transform to the bivariate sequences using Matlab’s discrete Fouriertransformbased“Hilbert”function. Forbothvaluesofcorrelationwecomputed ^ PLV sample , √ ^ PLV 2 ub sample , and ^ PLV circgauss as a function of sample size N. A total of 1,000 Monte Carlo trials were performed for each value of coherence and N, andsamplemeansandvariancescomputedforeachmeasureofPLV.Resultsareplotted inFig. 2.6. All three measures exhibit bias for small N. Note that while ^ PLV 2 ub sample is an unbiasedestimatorofPLV 2 ,itssquarerootisnotanunbiasedestimatorofPLValthough itdoesexhibitthelowestbiasamongthethreeestimators. Asexpected, ^ PLV ub sample also has the largest variance. It is also clear from Fig. 2.6 that the variance when esti- mating PLV from the coherence is significantly smaller than that of ^ PLV sample . Con- sequently, for data that are approximately Gaussian, using coherence to compute phase locking gives a more reliable (lower variance) interaction measure than does sample PLV. In comparison, in the case of the von Mises distribution, we saw that the sample PLVisamaximumlikelihoodestimatorwhichisasymptoticallyefficientandinpractice tendstoexhibitclosetominimumvariancebehaviorevenforsmallsamplesizes. 2.3.2 RoesslerOscillatorSimulations Roessler oscillators are commonly used models of weakly coupled stochastic oscilla- tors. We generated two Roessler oscillators 1 and 2 using the equations described in [SWD + 06]: 46 _ j = 0 B B B @ _ X j _ Y j _ Z j 1 C C C A = 0 B B B B @ ! j Y j Z j + [ ∑ i̸=j ϵ i;j (X i X j ) ] + j ! j X j +aY j b+(X j c)Z j 1 C C C C A (2.24) where i;j 2 f1;2g. The parameters were a = 0:5, b = 0:2, c = 10, ! 1 = 1:03, ! 2 = 1:01 and j is standard Gaussian noise. A range of values of were used as showninFig. 4.6. Parameterϵ i;j controlstheamountofcouplingfromthei th tothej th oscillator and is set such that ϵ 12 = ϵ 21 = ϵ implying a bidirectional coupling. In this study we consider two cases: coupling (ϵ > 0) and no coupling (ϵ = 0). For each of 1;000trialswegeneratedtimeseriesuptolength 9;000sampleswithasampleinterval ∆t = 0:02. WeusedthesimulatedtimeseriesX 1 andX 2 toinvestigatethebehaviorof ^ PLV sample and ^ PLV circgauss forthecaseswithandwithoutcoupling. Wegeneratedreceiveroper- ating characteristic (ROC) curves [Swe96] showing the fraction of true positive values (PLV> whenϵ> 0)versusfalsepositives(PLV> whenϵ = 0)asafunctionofthe threshold . Comparing the ROC curves for the two different PLV estimators we can determine which has the better sensitivity vs. specificity performance, as indicated by theareaunderthecurve,AUC(anAUC=1representserrorfreedetection). A sample ROC curve is shown in Fig.2.7(a) for the two estimators, indicating superior detection performance for ^ PLV circgauss compared to ^ PLV sample . Figs 7b-d 47 (a) (b) (c) (d) Figure 2.7: Roessler oscillator simulations (a) Comparison of ROC curves of ^ PLV circgauss (blue) and ^ PLV sample (red) whenϵ = 0:15 for coupled oscillators, stan- darddeviation = 1:5andnumberofsamplesL = 5000. (b)AreaunderROCcurveas a function of coupling parameter ϵ when = 1:5 and L = 5000. (c) Area under ROC curveasafunctionofnumberofsamplesLwhen = 1:5andϵ = 0:15. (d)Areaunder ROC curve as a function of the standard deviation of the noise when L = 5000 and ϵ = 0:15. show the AUCs for the ROC curves over range of (b) coupling parameters, (c) num- ber of samples, and (d) variances. Results consistently show superior performance for ^ PLV circgauss , particularly for small values of the coupling parameter and for a rela- tively small number of samples. As these parameters increase, results for both estima- tors are similar. Since the ^ PLV circgauss is a deterministic monotonic function of the 48 sample coherence, it follows that identical ROC curves would be obtained by replac- ing ^ PLV circgauss with the sample coherence. These results indicate that even in the cases where the data are not Gaussian, PLV calculation based on the Gaussian model (orequivalently,thesamplecoherence)cangivesuperiordetectionofcouplingthancan be achieved using the sample PLV. This seems to indicate that care should be taken in interpreting PLVs, since they do not necessarily provide different or more reliable insight into data than does the sample cross-correlation. This is the case even in this casewherethedataarenotobviouslyGaussian. Wenowturntothecaseforexperimen- taldata. 2.3.3 AnalysisofLFPdata Phase locking values have been used to analyze invasive (cortical and depth electrode) andnoninvasive(EEG,MEG)recordings. Inourworkweareinterestedinmodelingsig- nals ranging in scale from local field potentials from microelectrodes through invasive recordings with larger electrodes as well as EEG and MEG. Of these, microelectrode LFPshavethemostlocalizedsensitivity. Theotherrecordingscanbeviewedasequiva- lenttolinearmappingsofLFPsatthecorticallevel(orequivalently,thecombinationof dendritic and volume current sources that give rise to them) by integration with respect to the field sensitivities of the electrodes or magnetometers. These macro recordings should therefore be more Gaussian than the microrecordings through a law of large numbers argument. Conversely, multiunit and single unit recordings reflect the spiking activity associated with action potentials, and we would expect these to be highly non- Gaussian. Here we investigate the degree to which LFP recordings can be assumed to be Gaussian from the perspective of their relative phase distributions. In other words, weexplorewhetherPLVscomputedfromcoherenceproduceequivalentresultstothose computeddirectlyusingthesamplePLV. 49 For this study we used the LFP data described by [BCN93], which has been widely studied over the past two decades [Bre95, BDY99, DBYL00, BDL + 04]. Recordings were made using 51 m diameter bipolar electrodes separated by 2:5 mm in macaque monkeys. Measurementsfromatotalof15electrodepairsintherighthemisphere,with approximatelocationsshowninFig. 2.8,wereanalyzed. Intheexperimentthemonkey Figure 2.8: Locations of 15 electrode pairs in the right hemisphere (reproduced from [LDB01]) was trained to depress a lever and wait for a visual cue to either go (release) or no-go (not release). The go and no-go cues were a diamond and line pattern as shown in Fig. 2.9,withthechoiceofdiamondorlineasthegocuechangingbetweenexperiments. (a) (b) (c) (d) Figure 2.9: Four visual cues represented to the monkey in the experiment (a) right slantedline(b)leftslantedline(c)rightslanteddiamond(d)leftslanteddiamond. 50 Thecuewasgiven115msaftertheleverwasdepressed, andarewardgivenifthemon- key responded correctly to the go cue within 500ms. In the results presented below, we examine phase-locking at two time intervals: early response (120 25 msec after cue) andlateresponse(26025msecaftercue). Datawereoriginallysampledat200Hzand badtrialsdiscarded. Theremaining10;178trialscontained5225gotrials(gocueisline for2322trialsanddiamondfor2903trials)and4953no-gotrials(no-gocueisdiamond for2407trialsandlinefor2546trials). Forthisstudy,thedataforeachchannelandtrial wasfirst band-passfiltered to afrequencyband 13-30Hz byusing two-wayellipticIIR filtering of order 14. Hilbert transforms were then applied to each of these time series using Matlab’s “Hilbert” transform function. Because of the use of bipolar electrodes, the data should not be sensitive to volume conduction effects and therefore we restrict attentionheretoPLVanddonotcomputePLI. We first consider the relative phase distribution between pairs of electrodes. In Fig.2.10, we show goodness of fit between the empirical relative phase distribution extractedfromLFPdataandthetheoreticaldistributionsbasedonvonMisesandcircu- larlysymmetricGaussianmodels. Tocomputethetheoreticaldistributionsweassumed zeromeanandestimatedtheconcentration(forvonMises)andcoherence(forGaussian) parametersfromthesampledata. Theelectrodesandthetasksareselectedbasedonthe results represented in Table 2.2 which will be discussed later. The common feature of all selected electrodes is that they represent significant interaction when using sample PLV and circularly symmetric Gaussian PLV. We used a chi-square goodness of fit test for the two distributions [LS05]. We list the p-values for the selected pairs in Table 2.1. We were unable to reject the null hypothesis at = 0:05 for either distribution for anypairsshown,withtheexceptionthatthevonMiseswasrejectedbetweenelectrodes 3 and 5 at 120 msec as shown in Fig. 2.10(c) and listed in Table 2.1. We also plot 51 ElectrodePairs vonMises Gaussian 2and 3at 120msec 0.97 0.99 2and 5at 120msec 0.55 0.53 3and 5at 120msec 0 0.51 1and 7at 260msec 0.48 0.54 1and 8at 260msec 0.47 0.61 4and 7at 260msec 0.50 0.85 4and 8at 260msec 0.43 0.54 7and 8at 260msec 0.09 0.46 Table2.1: Thep-valuesofChi-squarestatisticforgoodnessoffitofthevonMisesand circularly symmetric Gaussian phase models. In only on case are we able to reject the nullhypothesis: electrodepair3and5forthevonMisesmodel. . ^ PLV sample versus ^ PLV circgauss using all pairwise combinations in Fig.2.11. We can concludefromtheseresultsthateitherdistributionisprobablyadequatetorepresentthe truephaserelationshipsinthis LFPdata. Wecomputed ^ PLV sample and ^ PLV circgauss for each electrode pair for each of the 18 experiments for the early and late intervals with trials sorted in two ways: (a) by visual cue (diamond vs. line) and (b) task (go vs no-go). Since the visual cues were switched in different experiments, the two different groupings allowus to differentiateinteractions associated with cue vs. those associated withtask. To assess significance of interactions we applied permutation testing by trial- shufflingwithinelectrodepairstoobtainanulldistribution(nointeraction). Usingthese distributions, computed separately for each pair and each condition, we then converted thePLVvaluestop-values. Weaccountformultiplehypothesestestingbythresholding withafalsediscoveryrate(FDR)of0.01. Foreachpairwethencomputedthenumberof experiments(outof18)inwhichPLVvaluesindicatedasignificantinteraction. Results are reported in Table 2.2 and illustrated in Fig. 2.12. Only pairs showing a minimum of5outof18experimentsineithercueortaskbasedcomparisonwithsignificantPLVs 52 areincludedinTable2.2whereasnothresholdwasappliedforrepresentationofresults inFig. 2.12. PLVtype EarlyResponse: LateResponse: sample Pair Diamond Line GO NO-GO Pair Diamond Line GO NO-GO 2-3 8 0 5 3 1-7 4 3 7 0 2-5 14 0 9 5 1-8 6 4 10 0 3-5 8 0 4 4 4-7 4 1 5 0 4-8 6 5 11 0 7-8 4 5 9 0 Gaussian Pair Diamond Line GO NO-GO Pair Diamond Line GO NO-GO 2-3 9 0 5 4 1-7 3 4 7 0 2-5 17 0 10 7 1-8 5 6 11 0 3-5 10 0 5 5 4-7 3 4 7 0 4-8 6 6 12 0 7-8 5 2 10 0 Table 2.2: Comparison of the number of experiments showing significant interactions betweenelectrodepairsforeachoffourconditions(Diamond,line,go,no-go)computed using sample PLV and Gaussian PLV for early and late response periods. Significance was assessed using FDR = 0.01 based on permutations of trial indices. Only pairs for which a minimum of 5 experiments show significant interactions in either cue or task basedcomparisonareshown. . In Table 2.2 we report the electrode pairs with significant PLVs for diamond vs. line and for go vs. no-go for both the early and late period. These results are included for twodifferentPLVestimators: ^ PLV sample ,whichistheMLestimatorforthevonMises distribution, and ^ PLV circgauss . Consider first the results for ^ PLV sample . For the early response we consistently see significant PLVs between electrode pairs 2-3, 2-5 and 3-5 for the diamond stimulus while there are none for the line. For this same early period, when sorting trials by go vs. no-go (both of which contain diamond and line stimuli) we see no indication that the go condition produces consistently more significant PLVs than the no-go condition. Note that electrode pairs 2-3, 2-5 and 3-5 are in striate/pre- striate cortex. If we now look at results for the late response we see that the above 53 observation is now largely reversed. In this case it is the go condition that leads to significantinteractionsbetweenstriate/prestriateandmotor/pre-motorcortices(pairs1- 7,1-8,4-7,4-8,7-8)whiletheno-goconditionshowsnosignificantinteractionsforany of these pairs. If the same data are sorted according to diamond vs. line we see no clear preference for significant interactions. These results might be expected since the earlyresponse(12025msecaftercue)shouldreflectvisualprocessingwhilethelater response (260 25 msec after cue) should also include visual/motor interactions. We now turn to the comparison in Table 2.2 between the two PLV estimators. While the numberofexperimentsshowingsignificantinteractionsisnotidenticalbetweenthetop and bottom halves of the table, the pairs of electrodes showing significant interactions are exactly the same for the two estimators, with very close agreement between the number of experiments with significant interactions for diamond vs. line in the early response and go vs. no-go for the late response. The Gaussian model consistently produces equal or slightly larger number of significant interactions in each column in thetable,possiblyreflectingsuperiordetectionpowerresultingfromthelowervariance oftheestimatorshowninFig. 2.6. 2.4 DiscussionandConclusion The primary purpose of this chapter was to explore the properties of the phase lock- ing value as it is frequently applied to electrophysiological data, i.e. as a measure of variabilityoftherelativephasebetweentwosignalscomputedbyaveragingovermulti- pletrialsandrelativelyshorttimewindows[MLDEE00,LRM + 00,VLRM01,BPFR01, DRKW08,OGT + 10]. ThesamplePLVisabiasedmeasure,aspreviouslyshownby[VVWW + 10]. While the degree of bias depends on distribution, the bias in the squared sample PLV was 54 shown to be a simple function of sample size N, equ. (2.10). Rearranging this result produces an unbiased estimator, equ. (2.11). It is interesting to note that this result is identical to Vinck’s PPC (Pairwise Phase Consistency) measure. In practice bias is relatively small for N > 50 so that the sample PLV should be adequate unless sample sizeissmall. The sample PLV is a maximum likelihood estimator when the data follow a von Misesdistributionimplyingefficiencyoftheestimator. WealsosawthatPLVisamono- tonic function of the concentration parameter and independent of the mean. Similarly we also saw that for jointly Gaussian data, the PLV is a monotonic function of coher- ence and also independent of mean. This latter result was based on marginalizing the jointdensityforthecomplexrepresentationofthetwosignalswithrespecttothesignal amplitudes leaving a distribution as a function of relative phase. In contrast to the von Mises case, the sample PLV is not a maximum likelihood estimator, and in fact we saw from Fig. 2.6 that computing PLV directly from the sample coherence leads to a lower varianceestimatethanthesamplePLV. It is not surprising that PLV depends only on coherence for Gaussian signals, since jointly Gaussian processes are completely characterized by their mean and covariance. Nevertheless, the consequence of this observation is perhaps less obvious: that if data are reasonably well modeled as jointly Gaussian, then the PLV provides no informa- tion that is not already contained in the coherence. We saw that this appears to be the case for the macaque LFP data analyzed above as well as the simulation based on the Roessleroscillator. ComparingsamplePLVresultswithaPLVcomputeddirectlyfrom the sample coherence we see very little difference between results in terms of the pairs of electrodes exhibiting significant interaction when using the two different measures. ThedistributionplotsinFig. 2.10,andassociatedhypothesistestsinTable2.1,indicate that both the von Mises and circularly symmetric Gaussian distributions are adequate 55 for this data. This supports the earlier statement that PLV does not add information not containedinthecoherence. It appears reasonable to expect that similar results to these LFP studies would be obtained when analyzing MEG and EEG data, since they share with the LFP a depen- denceprimarilyondendriticcurrentsandtheassociatedvolumeorconductioncurrents. Differences between these data depend on the lead field sensitivities of the transducers, with MEG and EEG being sensitive to far larger regions of cortex than the microelec- trodesusedtoacquirethedatainthestudy. Throughthelawoflargenumbers,wewould expect as we integrate over larger regions of cortical activation, that data would tend to becomemoreratherthanlessGaussian. It is important to note that the above observations are applicable only to data that can be considered approximately stationary over the interval over which the PLV is calculated. This is frequently the case in published studies, where either PLV is cal- culated over a short time window, or using a single time sample and averaging over trials[LCL03,LRM + 00,HRS04,RDK + 06]. IfPLViscomputedoverlongertimeperi- ods for which the behavior cannot be considered to be stationary, then the equivalence of PLV and coherence would not necessarily be retained. We also emphasize that the above arguments apply only to LFPs, EEG and MEG. Multiunit recordings are clearly non-Gaussianandwewouldnotexpecttoseeequivalenceinthiscase. In several instances false phase locking can arise as a result of linear mixing or cross-talk between time series: (a) use of common reference in microelectrode or EEG recordings, (b) EEG or MEG recordings with overlapping sensitivities in the lead field, and (c) cortical current density maps computed from EEG or MEG data in which the low resolution causes interference between sources. Since the PLV is independent of themeanphasedifference,thismeasurecannotdifferentiatebetweenazero-meanphase difference, which can be explained by linear mixing, and a non-zero phase lag, which 56 cannot be caused by linear mixing. The Phase Lag Index (PLI) was designed to be robust to this problem. In Figures 2.3 and 2.5 we show the dependence of PLI on the parameters of both the von Mises and Gaussian distributions. Unlike PLV, PLI is a function of the mean as well as the variance of relative phase. While this does produce the desired robustness to linear mixing, it is also a function of both mean and variance so that it is not possible to assess the strength of the interaction without also knowing the phase. One solution to this problem might be to compute both the mean and PLV of the relative phase, and use the mean to test for linear mixing and the PLV to assess the strength of the detected interaction. A similar problem is encountered when using coherence as an interaction measure: a real value of coherence can arise simply due to linear mixing. An effective solution to this problem is to use only the imaginary part of coherence [NBW + 04], which cannot be nonzero without true interaction. We note that even in the Gaussian case, where PLV and coherence are equivalent, the PLI and imaginary coherence are not: while both depend on the complex value rather than the magnitude of the coherence, they are different functions of this parameter. In the next chapter,wewillexploretheinteractionmeasuresinthepresenceoflinearmixing. 57 (a) (b) (c) (d) (e) (f) (g) (h) Figure 2.10: Goodness of fit between empirical distribution of relative phase extracted from LFP data and parametric von Mises and circularly symmetric Gaussian distribu- tions. Theconcentrationparameter(vonMises)andcross-correlationparameter(Gaus- sian) were directly estimated from the sample data. (a-c) pairs at 120 msec (d-h) pairs at260msec 58 Figure2.11: Scatterplotof ^ PLV vonmises versus ^ PLV circgauss computedfrommacaque LFPdata. PLVtype Earlyresponse: DiamondvsLine Lateresponse: GOvsNO-GO Sample Gaussian Figure 2.12: Phase synchronization networks constructed from the results in Table 2.2. Green: significant PLVs occur when diamond is presented. Magenta: significant PLVs occur when line is presented. Blue: significant PLVs occur during go condition. Red: significant PLVs occur during no-go condition. The thickness of the edges represents thenumberofexperimentssessionsinwhichsignificantsynchronizationoccurs. . 59 Chapter3 PartialPhaseSynchronization 3.1 Introduction In chapter 2, we explored the relationships between PLV and coherence. A significant limitation of PLV is that it is a bivariate measure which cannot differentiate between directandindirectinteractionsinamultiple-nodenetworksincetheinteractionsbetween neuronalpopulationsarenotpairwise,buttypicallyinvolvemultipleareasorelectrodes. Whencomputingpairwisecorrelations,interactionsbetweenpairscanbeproducedeven ifthe pairdoes nothaveadirect interaction butinstead bothof theareas areinteracting with a third. For this reason, multinode network models are sometimes inferred using a multivariate model in which partial correlations are computed to differentiate direct from indirect interactions between electrode pairs. These networks can be represented usingmultivariateGaussianmodelswherethepartialcorrelationsaregivenbythenon- zeroentriesoftheinverseofthecorrelationmatrix[Whi09,MMH00]. Recently Canolty et. al [CCK + 12] have developed a multivariate extension of the von Mises distribution that similarly is able to differentiate between direct and indi- rect interactions, but with respect to phase coupling rather than correlation. Schelter et al. [SWD + 06] investigated the inversion of the matrix of pairwise PLVs to compute a non-parametric estimate of partial PLV. Their approach is analogous to the inversion of cross-correlationandcross-spectralmatricestocomputepartialcorrelationsandpartial- coherence,respectively. 60 In this chapter, we derive an analytical expression for partial PLV under the mul- tivariate circular complex Gaussian model [Goo63]. We show that under this model, partial PLV is a nonlinear function of partial coherence, which is computed from both the phase and amplitude information of the signals. We demonstrate our partial PLV measure in detection of synchronization of Roessler oscillators. We further show that localfieldpotentialsfromamacaquemonkeystudyapproximatelyfollowourmodeling assumptions and our measure of partial PLV has reduced estimator variance compared tothenon-parametricPLVin[SWD + 06]. 3.2 Methods In this section, we discuss two different multivariate measures of phase synchrony. We first describe the non-parametric multivariate partial PLV measure proposed by Schel- ter et al. [SWD + 06]. We then derive a parametric measure of partial PLV using the multivariatecircularcomplexGaussianmodel. 3.2.1 Non-parametricPartialPLV Schelter et al. [SWD + 06] proposed inverting the matrix of pairwise PLVs to com- pute a non-parametric partial PLV. This is motivated by the analogous inversion of cross-correlation and cross-spectral matrices to compute partial correlation and partial coherence, respectively [Dah00]. First, the matrix of pairwise phase synchronization indices is constructed such that P mn = E [ e jϕm;n ] , where ϕ m;n ≜ arg ( zm(t)z n (t) jzm(t)zn(t)j ) , m;n2f1; ;Mg ,M is the total number of oscillators andfz m (t);z n (t)g is the pair ofanalyticsignalsofinterest. Then, inverting matrixP and normalizing with the diagonal elements leads to the partialPLV: 61 PPLV nonparam △ = (P 1 ) ab √ (P 1 ) aa (P 1 ) bb (3.1) between nodes a and b. Since the above expression is based on the sample bivariate PLV,itisanon-parametricestimateofpartialPLV. 3.2.2 Parametric partial PLV Using the Multivariate Circular ComplexGaussianModel In this section we introduce a parametric measure of partial PLV as an extension to the Gaussian based bivariate PLV given in section 2.2.3 using the multivariate circular complex Gaussian model. Given jointly Gaussian, zero mean, time series s i (t), i = f1; ;Mg,j = p 1,theanalyticrandomvectorz(t)isconstructedwithelements z i (t) =A i (t)e jϕ i (t) (3.2) where A i (t) = √ s 2 i (t)+ ~ s 2 i (t) and ϕ i (t) = tan 1 (~ s i (t)=s i (t)) are the instantaneous amplitude and phase of s i (t), respectively, and ~ s i (t) is the Hilbert transform of s i (t). Again we drop the time index t for convenience. The probability density function ofz iscircularcomplexGaussiandistribution[Goo63]: p(z) = 1 2 jK z j exp { z H K 1 z z } (3.3) where K 1 z is the inverse covariance matrix of z with mn-th entry mn e jmn , mn = nm > 0and < mn = nm <. For M > 2, letz = 2 4 x y 3 5 wherex = 2 4 z 1 z 2 3 5 andy = [ z 3 z M ] T ; then thedistributionofxconditionedony isalsocircularGaussianwithconditionalinverse 62 covariance matrixK 1 xjy = 2 4 11 12 e j 12 21 e j 21 22 3 5 as derived in App. B.1. Hence, the conditionaldistributioncanbewrittenas: p(xjy) =p(z 1 ;z 2 jz M ) = 1 2 jK xjy j exp [ x H K 1 xjy x ] (3.4) where M = f3; ;Mg. The rest is very similar to the theory in section 2.2.3. We first represent this distribution in terms of phases and amplitudes as p(ϕ 1 ;ϕ 2 ;A 1 ;A 2 jϕ M ;A M ). Further conditioning on A 1 and A 2 reveals the distribution ofthepair-wisephaseconditionedonalltheothervariables: p(ϕ 1;2 jA 1 ;A 2 ;A M ;ϕ M ) = 1 2I 0 (2 12 A 1 A 2 ) e 2 12 A 1 A 2 cos(ϕ 1;2 12 ) : As in the bivariate case, this is a von Mises distribution with coupling parameter 2 12 A 1 A 2 , verifying that the distribution of relative phase is a function of amplitudes and thus not independent of amplitude. As a final step we marginalize with respect to (A 1 ,A 2 )tocomputetheparametricpartialPLV: PPLV param △ = E [ e jϕ 1;2 jA M ;ϕ M ] (3.5) = p 2 ( 1 11 22 2 12 )[ w 3=2 2 F 1 ( 3 4 ; 5 4 ;1;w 2 ) + 3 4 w 5=2 2 F 1 ( 5 4 ; 7 4 ;2;w 2 )] (3.6) where w = 2 12 2 11 22 2 12 and 2 F 1 is the hypergeometric function. The derivation is the same as the derivation in bivariate case in equation 2.23. The data can be re-indexed in order to find the PPLV between any pair of nodes or channels. Since the signals have been band-pass filtered before computing their Hilbert transform, the inverse of theconditionalcovariancematrix,K xjy ,effectivelycorrespondstothepartialcoherence 63 matrix for that frequency indicating that the partial phase locking value for Gaussian time series is simply a function of the partial coherence, 12 = p 11 22 . To estimate PPLV param , we find the maximum likelihood estimate of the partial coherence and thanperformthetransformationin(3.6). Incontrast,wecomputecoherenceinorderto estimateGaussianbasedPLVwhichisimplicitlyabivariatemeasure. 3.3 Results 3.3.1 RoesslerOscillatorSimulations As in section 2.3.2 we used Roessler oscillators to test whether the aforementioned meauresofpartialPLVcanseparatedirectfromindirectphaserelationships. Wegener- ated 3Roessleroscillators i usingtheequationsasdescribedin[SWD + 06]: _ j = 0 B B B @ _ X j _ Y j _ Z j 1 C C C A = 0 B B B B @ ! j Y j Z j + [ ∑ i̸=j ϵ i;j (X i X j ) ] + j j ! j X j +aY j b+(X j c)Z j 1 C C C C A (3.7) wherei;j2f1;2;3g. Wesettheparametersa = 0:5,b = 0:2,c = 10,! 1 = 1:03,! 2 = 1:01, ! = 0:99 and j = 3:5 for all j and j is standard Gaussian noise. Parameters ϵ i;j controltheamountofcouplingfromthei th tothej th oscillatorandaresetsuchthat ϵ 12 = ϵ 21 = ϵ 13 = ϵ 31 = ϵimplyingabidirectionalcoupling. Wealsosetϵ 23 = ϵ 32 = 0, sothatthereisnodirectcouplingbetweenthe2ndand3rdoscillators. 64 Figure 3.1: (a) Estimated coupling between the 2nd and 3rd Roessler oscillators. For eachboxplot,thecentralmarkisthemedian,theedgesoftheboxarethe25-thand75-th percentiles,andthewhiskersextendtothemostextremedatapointsnotconsideredout- liers. The true coupling between is zero. (b) Q-Q plot of ranked Mahalanobis distance r 2 i versustheexpectedcorrespondingvalueof 2 3 fortheRoessleroscillators. We used the above equations to produce time series of 10;000 samples with sampling interval ∆t = 0:02, and then estimated the sample PLV (2.6), non-parametric par- tial PLV (3.1) and parametric partial PLV (3.6) between each pair of oscillators. By repeating the above procedure 1000 times, we computed the distribution of the above PLV measures. The corresponding box plots are shown in Fig. 3.1(a) for the coupling between 2 and 3 . Giventhatthetruecouplingiszerobetween 2 and 3 ,theparametric partialPLVisnotonlytheleastbiased,butalsohasthesmallestvariance. 65 (a) (b) Figure 3.2: ROC analysis of the Roessler oscillators (a) ROC curves for ϵ = 0:2 (b) AreaunderROCcurvesfordifferentvaluesofϵ(legendas(a)). To evaluate whether the signals X 1 , X 2 , and X 3 are jointly Gaussian, we used an extension of a univariate graphical procedure (Q-Q plot of a Gaussian distribution) to themultivariatecase[MM04]. TheMahalanobisdistance: r 2 i = (X i X) T S 1 (X i X) (3.8) whereX = [ X 1 X 2 X 3 ] T ,willbedistributedapproximatelyasachi-squaredistri- butionwith 3degressoffreedomiftheRoessleroscillatorsarejointlyGaussian. Multi- variatenormalityisrejectediftheQ-QplotoftheorderedMahalanobisdistancesversus the corresponding chi-square quantile is significantly non-linear. Fig.3.1(b) shows that the Roessler oscillators can deviate significantly from multivariate normality, depend- ing on the selected parameter ϵ. Despite the deviation from normality, the parametric partial PLV measure is able to detect the zero direct connection between 2nd and 3rd. The parametric partial PLV performs similarly to the non-parametric one when there is significantdeviationfromnormalityasshowninthesecondrowofFig.3.1. We also performed Receiver Operating Characteristic (ROC) analysis to compare the accuracy of methods in detecting coupled versus non-coupled nodes. For the same 66 Roessler system described above, we computed estimates of PLV between 2 and 3 (ϵ 23 = 0; non-coupled) and 1 and 2 or 1 and 3 (ϵ > 0; coupled). For each threshold value, we computed the false positive and true positive rates, producing the ROC curve given in Fig. 3.2(a). We then measured the area under the ROC curves for different valuesofϵ, as shownin Fig. 3.2(b). Both plots indicatethe non-parametric partial PLV hasbetterperformancethantheothertwomethods. 3.3.2 Localfieldpotentialsfromamacaquemonkey We performed a similar synchronization analysis to investigate oscillatory interactions of local field potential (LFP) time series, sampled at 200Hz, from a macaque mon- key implanted with transcortical bipolar electrodes at 15 sites in the right hemisphere (Fig. 2.8) while the monkey performed a GO, NO-GO visual pattern discrimination task [BCN93]. In this work, we used only 313 GO trials and focused on 120 msec and 265 msec after stimulus presentation in the frequency range of [10 15] Hz. Fig. 3.3 shows the Q-Q plot of ranked Mahalanobis distances versus 2 15 quantiles indicating that the data is approximately multivariate Gaussian at time instants 120 msec (early stimulus)and 265msec(responseonset). We then calculated channel-to-channel synchrony and repeated the procedure over 500 bootstrap samples of the data to estimate the variance of the PLVestimates. Figure 3.4 (a-c) shows the mean bootstrapped PLVs. The synchronization network from the bivariate sample PLV is quite different from the partial PLV ones. The synchronization networksfromthetwopartialPLVmeasuresaresurprisinglysimilar,howeverthestan- dard deviation of the parametric model is somewhat smaller than the non-parametric estimate(Fig. 3.4(d)). 67 (a) at 120msec (b) at 265msec Figure3.3: Q-QplotofrankedMahalanobisdistanceversusthecorrespondingquantile of 2 15 formacaqueLFPdata. (a) SamplePLV (b) ParametricpartialPLV (c) Non-parametricpartialPLV (d) Std nonparametric -Std parametric Figure 3.4: (a-c) Bootstrap mean of the PLV measures. (d) Difference of bootstrap standarddeviationofnon-parametricminusparametricpartialPLV. 68 3.4 Conclusion The use of multivariate phase synchrony measures is a valuable approach to detecting interactions in dynamic networks with multiple nodes. We have shown the qualita- tive difference between bivariate and multivariate approaches both in simulations with Roessler oscillators and macaque monkey LFP data. The parametric partial PLV mea- sure outperformed the non-parametric one in detecting non-coupled nodes when the data was approximately Gaussian (Fig. 3.1, top). Both measures performed similarly whenthedatawasnon-Gaussian(Fig. 3.1,bottom). ROCanalysisfurtherdemonstrated that the parametric partial PLV is better able to distinguish direct and indirect coupled Roessleroscillators. Since the parametric partial PLV is a non-linear function of partial coherence, does not constitute a novel interaction measure. Rather, it provides a new insight into phase synchronizationinmultivariatenetworks. WhenthetimeseriesareapproximatelyGaus- sian, which is often the case in experimental data, partial coherence and partial PLV containequivalentinformationaboutthesystemdynamics. 69 Chapter4 ConnectivityinthePresenceof Cross-talk 4.1 Introduction Asdiscussedinthepreviouschapters,coherenceandphaselockingvaluearecommonly used measures to investigate brain networks. We have shown that these two measures provide equivalent information when the signals are approximately Gaussian [APL13]. Therefore,inthischapterwewillusecoherenceasastandardinteractionmeasure. A challenging problem in functional connectivity analysis of EEG/MEG data is cross-talk or linear mixing of signals between measurement devices as a result of the diffusenatureoftheleadfieldsensitivityofthesemodalities. Neighboringelectrodesor sensors are sensitive to common sources within the brain resulting in spurious inter- actions between channels using standard interaction measures such as coherence or phase locking value [NSW + 97]. Mapping the EEG/MEG data into an estimated cor- tical source domain image [SP98] attempts to overcome this problem but a significant amountofcross-talkstillremainsbecauseofthelimitedspatialresolutionoftheresult- ing source maps. The point spread functions of these maps quantify the degree of mix- ing from each location to all other parts of cortex. Alternative methods such as beam- formers can also yield misleading results [SP98] although these can be ameliorated to some degree using the concept of a nulling beamformer [HPBL10]. This problem is termed “volume conduction” [NSS + 99] and have been increasingly addressed recently 70 [GHAEC08, GMDV + 11, HPBL10, MDGN08, NZN + 08, SG09, SND07]. Similarly, common reference electrode in EEG contributes a significant amount of interaction [FRBM88]. The problem of cross-talk generates false positive connections. Therefore, ithastobeovercomeforareliableinterpretationoffunctionalconnectivitynetworks. SeveralmeasuresincludingImaginaryCoherence(IC)[NBW + 04],PhaseLagIndex (PLI)[SND07]andLaggedCoherence(LC)[PM07]havebeenintroducedtoovercome this problem. These measures used the assumption of cross-talk being instantaneous or introducing negligible time lag due to quasi-stationary description of the forward model [SP98]. In other words, time lagged interaction is thought to be a better indica- tor of a true interaction than instantaneous interaction. It is well-known that the time lagged interaction causes an increase in the imaginary part of the complex coherence or the relative phase difference between signals. Based on this observation, Nolte et al. [NBW + 04] proposed to use IC ignoring the real part of the complex valued coherence whichcouldyieldspuriouslargevaluesduetocross-talk. ThePLIproposedbyStamet al. [SND07],ontheotherhand,quantifiestheamountofdeviationoftherelativephase aroundzero. Consistentrelativephasearoundzeroisrelatedtoinstantaneousinteraction thusanindicatorofacross-talk. Another measure proposed by Pascual-Marqui et al. [PM07] is the lagged coher- ence (LC). The interpretation of this measure is two-fold. First, it can be considered as a regression problem where the linear interaction between two complex processes is determinedbyalinearregressingcoefficient. Therealpartofthiscoefficientrepresents an instantaneous interaction whereas the imaginary part measures time-lagged interac- tion. LC quantifies the reduced variance of the residual of the linear estimation with an incorporation of an imaginary regressor in addition to a real regressor. In other words, an increased performance in the estimation with use of imaginary coefficient indicates 71 trueinteraction. Second,itisanimprovedversionofICwheretheeffectofcross-talkis subtractedinthedenominatorofIC. Recently, Hipp et al. [HHC + 12] has developed a new approach to identifying large scale interactions that overcomes the effect of linear mixing. The approach is based on anorthogonalizationprocedurethatremovestheinstantaneousinteractionbetweentime seriesandtheamplitudecorrelationisthenperfomedbetweentheseorthogonalsignals. Therefore,wecallthismeasure“orthogonalcoherence”(OC). Measures C, IC, PLI and LC are all bivariate approaches and do not consider the effect of the remaining signals in the multivariate network. Partial Coherence (PC) [Lau96,Whi09]isamultivariatemeasurewhichcomputesthecoherencebetweenresid- uals after removing the effect of other sources. This way, it quantifies whether or not two signals of interest are directly interacting with each other. Although PC does not necessarily consider the cross-talk artifact, it has also been used to reliably estimate functionalconnectivityinpresenceofcross-talk[MMH00]. Theprimarygoalofthischapteristoprovideacomprehensivereviewoftheexisting measureswhichareproposedtoovercomecross-talkproblem. WewillshowthatLCis invariant to cross-talk in a bivariate model of cross-talk which was indepedently shown in[EMZ + 12]andOCissimplyanotherinterpretationofLC.Inspiredbytheregressing on the other sources in PC, we formulated a new measure Partial Lagged Coherence (PLC) which exludes the cross-talk effect by using real regressors and computes LC betweentheresiduals. L1normregularizationisappliedontheregressioncoefficientsto avoid overfitting. We will show PLC is more robust against cross-talk in a multivariate framework where multiple cross-talks may occur. This approach is applicable to any electrophysiologicaldatawherecross-talkislikelytooccur. The organization of the chapter is as follows. In section 4.2, we first describe the cross-talkproblemandintroducethemathematicalnotation. Theconventionalmeasure 72 (a) (b) (c) Figure4.1: Differentinterferenceandlinearmixingschemes. Thesourcess 1 ands 2 are measured in channels x 1 and x 2 . (a) Ideal case: no linear mixing no interference. (b) linearmixingwithoutinterference. (c)linearmixingwithinterference. C is also recalled in this section. In the following section 4.3, we discuss two different modelsforcross-talk. Insection4.3.1,abivariatemodelisintroducedwhereonlysingle cross-talk is present. The measures IC, PLI, LC and OC are described in this section. Wethenformamultivariatecross-talkmodelwherePCisrecalledandPLCisproposed insection4.3.2. Intheresultssectionweexplorethemeasuresforbothmodels. Forthe multivariate model, we generate a realistic simulations using MEG model to evaluate themeasures. Wefinallydiscussourfindingsinsection4.5. 4.2 ProblemDefinition Supposeweareinterestedintheinteractionbetweentwocomplexvaluedsignalss 1 and s 2 on cortex as illustrated in Fig. 4.1. In practice, these signals are not observable but the recorded signals at the measurement devices are. Let’s denote these observed signalsasx 1 andx 2 asshowninFig. 4.1. Thecomplexrepresentationscanbeobtained either by short time Fourier Transform or by band-pass filtering followed by Hilbert transform. Inourwork,weusethelatterapproach. Intheidealcasewherethesystemis freeofcross-talkasshowninFig. 4.1(a),theinteractionbetweensourcess 1 ands 2 can 73 be directly estimated from the measurementsx 1 andx 2 using scale invariant measures. Coherence[NSW + 97]andphaselockingvalue[TRW + 98]arewidelyusedmeasuresto investigatelinearrelationshipbetweentwotimeseriesatagivenfrequency. Inourrecent work,wehaveshownthatthesetwomeasuresprovideequivalentinformationwhenthe signals are approximately Gaussian [APL13]. Therefore, the coherence can be used to estimate interaction between s 1 and s 2 given x 1 and x 2 . The coherence between zero meannarrow-bandcomplexsignalsisdefinedas C (x 1 ;x 2 )≜ jE[x 1 x 2 ]j √ E [ jx 1 j 2 ] E [ jx 2 j 2 ] (4.1) which is always less than or equal to 1 [NSW + 97]. We will call this quantity complex coherenceiftheabsolutevalueisnottaken. Thismeasurecanalsobeinterpretedasthe normalizedcross-correlationbetweentwonarrowbandcomplexsignals. Duetovolumeconduction,onesourcemightbepickedupatdifferentmeasurement devices; this is the problem of cross-talk (illustrated in Fig. 4.1(b)). Although source spaceanalysisprovidesmoreaccurateestimationofthesources,cross-talkproblemstill exists due to the fundamental ill-posedness of the EEG/MEG inverse problem. Either in sensor or source space analysis, this problem has to be overcome in order to reliably detectinteractionbetweensignals. In this study, we discuss the methods which are more robust to cross-talk. The problem can be extended to multiple signals with multiple cross-talks as represented in Fig. 4.1(c). We assume that cross-talk due to volume conduction does not introduce time shifts in the signals as also assumed by Nolte et al [NBW + 04]. This assumption stems from the quasi-static approximation of the forward problem. Therefore, all the linear mixing coefficients fa mn g in figures 4.1(b) and 4.1(c) are assumed to be real. 74 Note that the signalsfx 1 ;x 2 ;x 3 g represent the mixed signals either in sensor or source domain. In the next subsections we will describe different approaches proposed to reliably estimatebrainnetworksandweproposeanovelmeasurewhichislessaffectedbymul- tiplecross-talks. 4.3 Methods 4.3.1 InteractionMeasuresinPresenceofaSingleCross-talk A cross-talk model between two complex valued sources s 1 and s 2 is depicted in Fig. 4.1(b). Mathematically,wehave: x 1 (t) = s 1 (t)+a 12 s 2 (t) (4.2) x 2 (t) = a 21 s 1 (t)+s 2 (t) (4.3) where the inequalities a 12 < 1 and a 21 < 1 hold for real cross-talk coefficients so that the largest amount of power of the k-th source remains in the k-th recording sensor. The amount of cross-talk is determined by these cross-talk coefficients a 12 and a 21 . In thiscase,thecoherencewillyieldaspuriouslyincreasedvalue. Inordertodescribethe effects of cross-talk on coherence, we write the numerator of the coherence in equation 4.1(whichissimplythecross-correlation)betweenx 1 andx 2 intermsofs 1 ands 2 as: E[x 1 x 2 ] = E[(s 1 +a 12 s 2 )(a 21 s 1 +s 2 ) ] (4.4) = a 21 E [ js 1 j 2 ] +a 12 E [ js 2 j 2 ] | {z } Real +E[s 1 s 2 ]+a 12 a 21 E[s 1 s 2 ] | {z } Complex : 75 It can be seen that only the complex part is affected by true interaction. In other words, therealpartofthecross-correlationmaybelargeevenfornon-interactingsources. The imaginarypartofthecomplextermswillvanishifthereisnophasedifferencebetween sourcesandthephasedifferenceimpliestime-laggedinteraction. Basedonthisobserva- tion, several interaction measures have been proposed. We will discuss these measures inthefollowingsubsections. ImaginaryCoherence Sincetheimaginarypartofthecomplexcorrelationorthecoherencebetweenmeasure- mentscanbealteredonlybytrueinteraction,Nolteetal. [NBW + 04]proposedIC: IC (x 1 ;x 2 )≜ jIfE[x 1 x 2 ]gj √ E [ jx 1 j 2 ] E [ jx 2 j 2 ] : (4.5) where I denotes the imaginary part. This is simply the imaginary part of the com- plex coherence. The fundamental reasoning of the use imaginary coherence is that its nonzerovaluecannotbecausedbyvolumeconduction. PhaseLagIndex In order to measure consistency of the non-zero valued phase difference as a conse- quence of the time-lagged interaction, Stam et al. [SND07] has introduced PLI which quantifiestheasymmetryofthedistributionoftherelativephasearoundzero: PLI(x 1 ;x 2 )≜jE[sign(ϕ 1 ϕ 2 )]j (4.6) 76 where ϕ 1 and ϕ 2 are the phases of the complex signals x 1 and x 2 , respectively. The PLI ranges between 0 and 1. The fundamental motivation here is to discard phase dif- ferences around 0mod. Large values of PLI imply presence of consistent non-zero phasedifferencewhichcannotbeexplainedbyvolumeconduction. Whilethelowval- ues imply either no phase coupling or existence of volume conduction. The standard phasesynchronizationmeasuresofphasesynchronizationyieldlargevaluesinthelatter case. ThismakesPLImorerobustagainstvolumeconduction. LaggedCoherence Pascual-Marquietal. [PM07]hasintroducedtheLC,whichisdefinedas: LC (x 1 ;x 2 )≜ IfE[x 1 x 2 ]g √ E [ jx 1 j 2 ] E [ jx 2 j 2 ] (RfE[x 1 x 2 ]g) 2 : (4.7) ThisisanimprovedversionofIC.Herethenormalizationismodifiedtoaccountforthe cross-talk component of the similarity between the two signals. As shown in App. C.1 andindependentlybyEwaldetal. [EMZ + 12],thismakesLCindependentofthedegree of cross-talk and solely a function of the true interaction between the two signals. In otherwords,LCbetweenx 1 andx 2 isidenticaltoLCbetweens 1 nds 2 . Thisprovesthe cross-talk invariant property of the LC which is the fundamental reason of our interest onthismeasure. In order to better understand LC, we visit the concept of the best linear predictors. Suppose the samples of the zero mean complex random variablex 1 are estimated from thesamplesofanotherzeromeancomplexrandomvariablex 2 suchthat x 1 (t) =ax 2 (t)+e(t) (4.8) 77 where a is a complex linear coefficient which determines the linear relation between x 1 and x 2 and e(t) is the error term in the estimation. The real part of this coefficient can account for the amount of signal x 2 (t) that is linearly mixed into x 1 (t) (with zero phase lag). Therefore, by constraininga to be real one can estimate the cross-talk from x 2 tox 1 . In contrast, the imaginary part of this coefficient is related to the time-lagged interaction. ThesquaredLCcanthenbedefinedas 2 LC ≜ 1 min a2C E [ jx 1 ax 2 j 2 ] min a2R E [ jx 1 ax 2 j 2 ] (4.9) where the numerator of the quotient is the residual variance of x 1 predicted from x 2 using a complex coefficient and the denominator is the residual variance using a real coefficient. In this way, LC measures whether the residual variance is reduced with the incorporation of the time-lagged interaction, removing any effects of instantaneous interaction. Ifso,thisimpliescouplingthatcannotbecausedbycross-talk. Usingdifferentconstraintsforthespacesoftheregressioncoefficientinthenumera- toranddenominatorofthequotientmightbecomputationallydifficult. Toreducecom- plexity,wecanwriteeachcomplextermintermsofrealandimaginarycomponents: x R 1 (t)+jx I 1 (t) = ( a R +ja I )( x R 2 (t)+jx I 2 (t) ) (4.10) + ( e R (t)+je I (t) ) (4.11) = ( a R x R 2 (t)a I x I 2 (t)+e R (t) ) (4.12) + j ( a R x I 2 (t)+a I x R 2 (t)+e I (t) ) (4.13) 78 where superscripts R and I denote real and imaginary terms, respectively and i is an abbreviationfor p 1. Equatingtherealandimaginarypartsyields 2 4 x R 1 (t) x I 1 (t) 3 5 = 2 4 x R 2 (t) x I 2 (t) 3 5 a R + 2 4 x I 2 (t) x R 2 (t) 3 5 a I + 2 4 e R (t) e I (t) 3 5 : (4.14) Wedefine ~ y(t)≜ 2 4 y R (t) y I (t) 3 5 and y(t)≜ 2 4 y I (t) y R (t) 3 5 foranycomplextimeseriesy(t). Henceequation4.14canbewrittenas: ~ x 1 (t) = ~ x 2 (t)a R + x 2 (t)a I +e(t): (4.15) LC quantifies the amount of improvement in the estimation by involving the term with a I inadditiontotheonewitha R . Thesquaredlaggedcoherenceisthendefinedas 2 LC (x 1 ;x 2 )≜ 1 min a R ;a I 2R E [ ~ x 1 ~ x 2 a R x 2 a I 2 2 ] min a R 2R E [ ∥~ x 1 ~ x 2 a R ∥ 2 2 ] (4.16) where the numerator of the quotient is the variance of x 1 predicted from x 2 using a complex coefficient and the denominator is the variance using a real coefficient as in equation4.9butherewerepresentallthevariablesinrealdomain. OrthogonalCoherence Hippetal. [HHC + 12]hasproposedanorthogonalizationprocedurethatdiscountsspu- rios correlation caused by volume conduction. For simplicity, we call this measure 79 “orthogonal coherence” (OC). Mathematically, this method computes the coherence betweenthe“orthogonal”signalsas: OC (x 1 ;x 2 )≜ C (x 1 ;x 2?1 ) (4.17) wherex 2?1 =x 2 RfE[x 1 x 2 ]g E[jx 1 j 2 ] x 1 andRdenotestherealpartsothattherealpartofthe coherence between x 1 and x 2?1 is zero. Note that the term “orthogonal” is not used in itsgeneralmeaningbutthezeroinstantanouscorrelationbetweensignals. Weprovethat the OC is equivalent to the LC in App. C.2. Therefore, we will not include the results ofthisapproachinthispaper. 4.3.2 InteractionmeasuresinPresenceofmultipleCross-talk In contrast to the single cross-talk model shown in Fig. 4.1(b), there may be multiple sources causing cross-talk as depicted in Fig. 4.1(c). This model better mimics the behaviorofelectrophysiologicalsignals. Mathematically,wehave: x 1 (t) = s 1 (t)+a 12 s 2 (t)+a 13 s 3 (t) (4.18) x 2 (t) = a 21 s 1 (t)+s 2 (t)+a 23 s 3 (t) x 3 (t) = a 31 s 1 (t)+a 32 s 2 (t)+s 3 (t) where the linear mixing coefficients are real so that the cross-talk is instantaneous and fa 12 ;a 13 ;a 21 ;a 23 ;a 31 ;a 32 g< 1sincek-thsensorisspatiallyselectivetothek-thsource fork2f1;2;3g. Thismodelcanbeextendedtothecaseofmorethanthreesources x(t) =As(t) (4.19) 80 where the linear mixing matrix A is presumably a diagonally dominant and sparse matrix. The sparsity assumption follows from the fact that most of the sources have negligible to no effect on cross-talk because the effect drops off as the square of the distance from the two nodes of interest [NSW + 97]. In the following, we consider a well-knownmultivariatemeasurePartialCoherence(PC)andPartialLaggedCoherence asanovelmodificationofLC. Partialcoherence Partial coherence or correlation (PC) is a measure of direct interaction between two processes after excluding the effects of the set of remaining processes [Lau96, Whi09] andhasbeenusedby[MMH00]todimishthecross-talkartifact. PCcanbeviewedasa correlationbetweentworesidualsafterregressingoutthesignalsfromothernodes. This methodallowsustoobtainameasureofdirectinteractionbetweentwonodeswhichcan notbeexplainedbythelinearcontributionsfromothernodesinthenetwork. Theseother sourcescausingcross-talktothetwosourcesofinterestarecalledinterferers. Say the interferers are denoted by a vectory where each entry represents a signal causingcross-talkandweregressthesignalsofinterest,say,x k andx m ony toremove thelineareffectfromtheinterferers: ^ x k (t) = c T k y (4.20) ^ x m (t) = c T m y (4.21) 81 wherethecomplexvaluedvectorsfc k ;c m g2C N2 representsetsofregressioncoeffi- cientsforeachequationandN isthetotalnumberofnodesinthenetwork. Theresiduals determinethelinearcross-talkeffectsoftheinterferersonthesignalsofinterest: r k (t) ≜ x k (t) ^ x k (t) (4.22) r m (t) ≜ x m (t) ^ x m (t): (4.23) Thesecomplexvaluedvectorsfc k ;c m gareobtainedbyminimizingthevariancesofthe residuals ^ c k = argmin c k 2C N2 E [ jr k j 2 ] = argmin c k 2C N2 E [ x k c T k y 2 ] (4.24) ^ c m = argmin cm2C N2 E [ jr m j 2 ] = argmin cm2C N2 E [ x m c T m y 2 ] : Thepartialcoherencebetweenx 1 andx 2 afterremovingtheeffectsofinterferersisthen definedas: PC ( x k ;x m ;fx i g i̸=k;m ) ≜ jE[r k r m ]j √ E[jr k j 2 ]E[jr m j 2 ] (4.25) which has limits between 0 and 1. PC can readily be calculated from the inverse of the coherencematrix,alsocalledconcentration,orprecision,matrix[Whi09]. However,the matrix inversion of ill-conditioned coherence matrices may lead to incosistency prob- lems. Numericalstableestimationofthecoherencematrixisobtainedbyregularization approaches [BGdN06, FHT08]. Here we describe the usage of regularization approach in terms of regression terms. We apply regularization method least absolute shrinkage and automatic variable selection (lasso) on the regression coefficientsfc k ;c m g so that simultaneousshrinkageandautomaticvariableselectionareperformed[Tib96]. Thisis done by forcing L1 penalty which consistent with the sparsity assumption of the linear 82 mixingmatrixAinequation4.19. TheestimationoftheregressioncoefficientswithL1 penaltyisthendescribedas ^ c k = argmin c k 2C N2 { E [ x k c T k y 2 ] + k ∥c k ∥ 1 } (4.26) ^ c m = argmin cm2C N2 { E [ x m c T m y 2 ] + m ∥c m ∥ 1 } where k and m aretuningparameterscontrollingthelevelofsparsityinf^ c k ;^ c m g. The residualsandPCarecomputedinthesamewayasinequations4.23and4.25. PartialLaggedCoherence WeformamultivariateextensionofLCbyremovingthecrosstalkeffectscausedbythe interferingsignals. ThismethodissimilartothePCexcepttheemphasisontimelagged interaction. We remove only the instantaneous, i.e. cross-talk, linear effects from the interferingsignalsandcomputetheLCinsteadofCbetweentheresidualsofthesignals of interest. To do this, we assume the regression vectorsfc k ;c m g in equation 4.21 are real valued so only the instantaneous effects will be obtained and we focus on the time lagged interaction by using LC. In practice, we do not expect cross-talk effect occur at everychannel. Therefore,weencouragesparsityconstraintonregressioncoefficientsby applyingL1normpenaltysimilartotheregularizationinPC.Theregressioncoefficients arecalculatedbysolvingtheequations: ^ c k = argmin c k 2R N2 { E [ x k c T k y 2 ] + k ∥c k ∥ 1 } (4.27) ^ c m = argmin cm2R N2 { E [ x m c T m y 2 ] + m ∥c m ∥ 1 } 83 where k and m are tuning parameters as in equation 4.26. Partial lagged coherence (PLC)isthencomputedbymeasuringtheLCbetweenresiduals: 2 PLC ( x k ;x m ;fx i g i̸=k;m ) ≜ LC(r k ;r m ) (4.28) = 1 min a R ;a I 2R E [ ~ r k ~ r m a R r m a I 2 ] min a R 2R E [ ∥~ r k ~ r m a R ∥ 2 ] : wherethevectors~ r fk;mg and r fk;mg calculatedfromr fk;mg aredescribedinsection4.3.1. SelectionofTuningParameters Tuning parameters in equations 4.26 and 4.27 control the amount of regularization. Therefore, good choice of these parameters is crucial for an accurate linear estimation. We use k-fold cross-validation to estimate prediction error for different regularization parameters. The data is partitioned into k “folds” where k-1 folds are used as training data and and the rest one fold is used for validation. The process is repeated k times and mean and standard deviation of each model are estimated. Hastie et al. [HTF + 09] described“onestandarderrorrule”toselectthemodelwhichisthelargestvaluesuch thattheerroriswithin1standarderrorofminimummeansquarederror. Thismethodis implementedinaMATLABpackagecalled lasso.m. 4.4 Results We used two different cross-talk models to analyze the methods described above. In the first model, only a single cross-talk is studied which was illustrated in Fig.4.1(b). Bivariate approaches such as C, IC, PLI and LC are compared for this model. More realistic model with multiple cross-talks as shown in Fig. 4.1(c) is then studied where PCandPLCarealsoincludedintheanalysis. 84 4.4.1 SingleCross-talk Inordertostudythecrosstalkeffectontheinteractionmeasures,wesimulatedabivari- ate network distorted by a cross talk. We generated bivariate zero mean Gaussian time seriesatfrequencyrangeof13-30Hz. Eachtimeseriesis5secondslongwithsampling frequency 200 Hz. Matlab’s Hilbert function is used to obtain complex representation of the signals. The coherence between the unmixed signals (represented ass 1 ands 2 in equation 4.3 ) is set to 0:2+j0:2 so that the time lagged interaction between signals is satisfied. The observed signals (represented asx 1 andx 2 in equation 4.3) are generated byrandomlymixingthesourcesbyvaryingrealcoefficientsa 12 anda 21 asdescribedin equation 4.3. We compute each interaction measure between mixed signals x 1 and x 2 andplotitasafunctionofcross-talkcoefficients. TheresultsshowninFig. 4.2indicate that all the measures except the LC are affected by the amount of cross-talk. In other words,theLCisinvarianttothecrosstalkcoefficientsasconsistentwiththetheoretical resultsprovidedinApp. C.1. 4.4.2 MultipleCross-talk Inpractice,morethanasinglecross-talkisthoughttobepresentinelectrophysiological recordings. Thus, we performed realistic simulations where the ground truth is known. We simulated multivariate time series with the same specifications as described in the previoussection4.4.1. Hilberttransformisagainusedtoobtaincomplexrepresentation of the signals. The interactions between the time series are predefined and the pres- ence of time lagged interactions is satisfied via non-zero valued imaginary coherence. These time series are used to model 8 focal dipole sources and their extended patches (approximately 5 cm2) on cortex. We modeled the geometry of Elekta Neuromag 306 85 (a) (b) (c) (d) Figure4.2: Interactionmeasuresasafunctionofcross-talkcoefficientsinpresenceofa singlecross-talk. (a)Coherence(b)Imaginarycoherence(c)Phaselagindex(d)Lagged coherence. Except from lagged coherence (d), all other interaction measures have high dependenceontheamountofcross-talk. channelsystemwhichconsistsof204planargradiometersand102magnetometers. The forward model is estimated via boundary element methods using openMEEG package [GPO + 10]. MEG electrodes are simulated by applying forward model to the simu- lated sources followed by addition of white Gaussian noise with a covariance matrix generated from real MEG sensor data. The SNR is defined as the ratio of a power of the signal in sensor domain to that of the total noise. We used minimum norm recon- struction method [MLL92] using BrainStorm [TBM + 11] for the inverse model based on overlapping spheres head model [HML99]. The estimated sources are computed by applying the inverse model to the simulated data in sensor space. By using different 86 Figure 4.3: Eight extended sources are selected on the cortical surface based on the locationsofimplantedelectrodesinamacaquemonkeyexperiment. assumptions on the head model for forward and inverse models, we avoid excessively optimistic expectations about the performance of the methods which is also called an “inversecrime”. Fig. 4.3showsthelocationsofthesimulatedextendedsourcesoncortex. Theassigment oftheselocationsisdeterminedbyanapproximatemappingfromtheoriginalplacement of the electrodes on the cortical surface of a monkey to a human cortical surface as describedin[APL13]. We first evaluate the methods on a long range interaction. We generated the simu- lations so that only the regions 2 and 7 (in Fig. 4.3) have non-zero time-lagged cross- correlation (coherence 0:2 +j0:2) and the rest of the pairs have zero cross-correlation. The SNR for each time series is set to 40. We computed each interaction measure 87 between estimated sources and repated the process 100 times. Fig. 4.4 shows the aver- age values connectivity matrices for each method over 100 Monte Carlo simulations. Sinceallthemeasuresaresymmetric,onlytheupperhalvesoftheconnectivitymatrices areshown. Accordingtothefigure, allmeasuresyieldseveralspuriousconnectionsbut thesevaluesareminimumforPLC. HistogramEstimation We measure the performance of each method in distinguishing connected pairs from unconnectedpairsbyanalyzingthedistancebetweenthedistributionsofeachtype. The distributionoftheconnectedpair(whichistheregions2and7inthiscase)isestimated over 100 Monte Carlo simulations whereas the distributions of the rest 27 unconnected pairs (all the pairs except the regions 2 and 7) is treated as a single distribution and estimated over 100 Monte Carlo simulations. In other words, we used 100 observa- tions for connected edges and 2700 observations for unconnected edges. We compute the normalized histograms via dividing the count of elements in each bin by the total number of bins. Fig. 4.5 shows normalized histograms of both connected edges and unconnected edges for each method. Larger distance between the histograms indicates a better performance in distinguishing connected edges from unconnected edges. The standard deviation of the connected edges in C and PC are small compared to that of unconnected edges. This will result in large number of spurious connections. Large value of a standard deviation of the unconnected edges can be explained by the large valuesofthesemesurescantake. Theseparationbetweendistributionsinthesecasesis notveryclearhencewillresultaweakperformance. IC,PLIandLCmayyieldabetter 88 (a) (b) (c) (d) (e) (f) Figure 4.4: Average values of measures over 100 Monte Carlo simulations where only the sources 2 and 7 are connected. Since the measures are symmetric only the upper halves of the connectivity matrices are shown. (a) Coherence (b) Imaginary Coherence (c) Partial Coherence (d) Phase Lag Index (e) Lagged Coherence (f) Partial Lagged Coherence 89 (a) (b) (c) (d) (e) (f) Figure4.5: NormalizedHistogramsofconnectedandunconnectededges performancethanCandPCastheprobabilityofhighvaluesfortheunconnectededges is mostly smaller than that of small values. The distinction between two types is more obvious in the case of PLC. The probability of high values for the connected edges and 90 the probability of low values for the unconnected edges are larger which is the desired conditiontodistinguishtwotypesofconnections. ReceiverOperatingCharacteristicCurves The performance of each method is illustrated by the behavior of receiver operating characteristic (ROC) curves [Swe96] showing the fraction of true positive connections (strength of interaction measure greater than when there is a true connection) versus false positive connections (strength of interaction measure greater than when there is no connection) as a function of the threshold . For instance, for a given value of , true positive rate is the ratio of the number of elements where the connection between regions2and7isgreaterthan to100andthefalsepositiverateistheratioofthenum- ber of elements where the connection between the rest of the pairs which are greater than to 2700. Comparison of ROC curves shown in Fig. 4.6 allow us to determine whichmethodhasthebettersensitivityversusspecificityperformance. Forareasonable amountoffalsepositiverate(FPR< 0:1),PLCyieldsthehighestvaluesofTPRindicat- ingabetterperformance. PLIperformsbarelybetterthanarandommethodandLCand ICperformverysimilartoeachother. CandPC,ontheptherhand,areunabletodetect truepositiveconnectionswhentheFPRissmallaswereobviousfromthehistogramsin Fig. 4.5. MultipleConnectionsandFractionalAreaUnderROC The area under ROC curve reflects the ability of a measure to distinguish a connected pair from an unconnected pair. The performance of a measure is thought to be better as this value approximates to one. In this work, we compute the area under ROC for the values where FPR is smaller than or equal to 0.1. Since this value is equal to 0.1 91 Figure4.6: ComparisonofROCcurvesofdifferentinteractionmeasuresinpresenceof asingleconnectioninamultivariatenetwork. in the ideal case, we normalize the values by 0.1 and call this measure “fractional area under ROC curve” (FAUC). We compute FAUC of each method for different schemes of connectivity networks as shown in Fig. 4.7. Fig. 4.8 shows FAUC values for these randomnetworks. PLCyieldsthe highestFAUCexceptforthenetworkwith6connec- tions. PC yields values very close to zero for the networks with connections 6, 10 and 12. Furthermore,theperformancesofCandPLIarenotasgoodasthoseofLCandIC. 92 (a) (b) (c) (d) (e) (f) Figure4.7: Truenetworksbasedonarandomselection 4.5 DiscussionandConclusion In this chapter we proposed a new method called PLC for the functional connectivity analysisinpresenceofmultiplecross-talks. TheLCistheoreticallyandviasimulations 93 Figure4.8: FractionalareaunderROCasafunctionofnumberofconnections showntobeinvarianttothesinglecross-talkasopposedtotheothermeasuresC,ICand PLIasdemonstratedinFig. 4.2. However,theseapproachesdonotaccountfortheaddi- tional cross-talk artifacts caused by other signals from different nodes in the network. PLCovercomesthecross-talkproblembyregressingouttheinstantaneouscontribution from interferers and focusing on the time lagged interaction between remaining resid- uals of the signals of interest. The focus on the time-lagged interaction is common in measures IC, PLI, LC and PLC since the cross-talk is thought to be instantenous. PLC differsfromLCintheregressionontheinterfererswhichisnotcountedinLC.Applying all these approaches to realistically simulated MEG data, we observed that PLC yield better results than the other approaches. MEG model is used to mimic the cross-talk problem but the same principles hold for EEG model as well. Furthermore, these find- ingscanbeappliedtomicroelectroderecordingswherecross-talkmightoccurbetween nearbyelectrodesorasaresultofacommonreferenceelectrode. 94 Identificationofconnectivitycanbeviewedasadetectionproblemwheretwotypes of errors could occur. (i) Type I errors where the interaction is falsely accepted when there is no interaction. (ii) Type II errors where the interaction is falsely rejected when there is an interaction. Based on these two types of errors we computed false positive and true positive rates. Although the measures discussed here attempt to eliminate the cross-talk, it can not be completely removed and thus the presence of false positive connectionsareinevitable. ThiswasobviousinFig. 4.4. Forinstance,inthesinglelong range connection in a multivariatemodel, none of the abovemodels were able to truely rejectalltheunconnectededgeswhiletruelyacceptingthetrueconnections. Therefore, we interpret the less number of false positive connections as a better performance for the same true positive rate. In Fig. 4.5, we plot the histograms of connected edges and unconnectededgeswherethelatteriscomputedovertherest27unconnectededges. We then used ROC curves to interpret these results. PLC demonstrated a better detection ratethantheotherapproachesascanbeseeninFig. 4.6. Reliable identification of a single long range connection is valuable but more than a single connections may be present in practice. Therefore, we randomly generated networks (shown in Fig. 4.7) where multiple connections were present. A fractional areaunderROCcurve(4.8)isusedfortheevaluationofthemethods. PLCdemonstrated abetterperformancethantheotherapproachesinmostofthecasesandLCandICyield very similar results. This indicates that there is not a significant difference between IC and LC in presence of multiple cross-talks. On the other hand, C, PC and PLI perform verypoorlyinmostofthecases. All the measures except the PLI quantifies the linear relationship between signals. These measures may not be able to capture the underlying network if the system is nonlinear. A construction of a nonlinear measure robust to cross-talk is beyond the scope of this work. We should note the limitation of PLC which is the requirement of 95 the definition of region of interest where the interactions and cross-talks are likely to occur. L1 norm penalty will reduce the possibility of the problem of overfitting in the regression analysis in case of large number of regions. Altough PLC is an extension of LC, it is inspired by PC in regressing out the interfering signals. Unlike PC, PLC can not be considered as a direct measure as it only removes the cross-talk effect from the other signals instead of the entire linear effect. In other words, PLC will not be able to distinguish two directly interacting sources from indirectly ones in which case PC will also fail in presence of cross-talk. In the analysis of LC, IC, PLI and PLC, it is possible that we may not capture the true instantenous interaction. In fact, we trade off the interaction with a negligible time lag for a larger time lag. This will cause an increaseinthefalsenegativeratewhilecontrollingthefalsepositiverate. The problem of cross-talk makes the brain connectivity analysis of EEG/MEG data verychallenging. Acarehastobetakenintheinterpretationoftheconnectivityresults. This indicates the value of studies trying to overcome this problem. We demonstrated that even in the simple model where only two signals are present the measures except LC may give spurious interaction. Therefore, we believe that the PLC approach we introducedhereaddressestheproblemdespiteitslimitations. In conclusion, we provided a comprehensive description of the methods proposed toestimatefrunctionalconnectivityinpresenceofacross-talk. Amongthesemeasures, LC was the most robust metric against cross-talk in the bivariate model. The measures C, IC and PLI are shown to be dependent on the amount of cross-talk. Our proposed approach PLC demonstrated additional robustness in presence of multiple cross-talks in a multivariate network. PLC combines LC with regressing interfering signals using realregressioncoefficients. Hence,itmeasuresthelaggedcoherenceafterexcludingthe cross-talkeffectfromothersources. Webelieveourapproachwillcontributetoreliable analysisoffunctionalconnectivityfromEEG/MEGdata. 96 Chapter5 ConclusionandFutureDirections 5.1 Conclusion In this dissertation we explored the functional connectivity analysis of the brain from electrophysiological signals. We presented solutions to three problems in the field: (1) thestatisticalrelationshipbetweenphasesynchronizationmeasures(2)theextensionof phasesynchronizationmeasurestothemultivariatecase(3)improvementofconnectiv- ityanalysisinpresenceofcross-talkinEEG/MEG. ThebivariatemeasuresPLVandcoherencearetwoofthemostcommonlyusedmea- sures to quantify functional connectivity. While PLV uses only the phase information of the signals, coherence incorporates amplitude information in addition to the phase information. Inchapter2wepresentedadetailedstatisticalanalysisofPLV.Wederived a closed form expression for this measure when the signals are Gaussian and proved that it is simply a monotonic function of another well-known measure coherence. The relationsofPLVtotheotherrelativephasebasedmeasuressuchasPLIandPPCarealso analyzed. Weconcludedthatphasesynchronizationmeasurescanbesummarizedintwo categories: (1)coherencebased(2)samplePLVbased. Thetheoreticalfindingsshowed one-to-one relationship between these two categories when the signals are Gaussian. We verified these findings via local field potential recordings from a macaque-monkey experiment. We believe that this conclusion will somehow reduce the confusion of the useoftherichnumberofconnectivitymeasuresinfunctionalconnectivitystudies. 97 In a multivariate network, bivariate measures such as coherence and PLV can not distinguish a direct interaction from an indirect interaction. Multivariate measures includingpartialcoherenceorpartialphaselockingvalueareproposedtocapturedirect interaction between two signals of interest. They attempt to understand if two signals are directly interacting or indirectly interacting through another signal. In chapter 3 we extended the bivariate PLV for Gaussian signals to the multivariate case and com- pared to the non-parametric approach. We compared the results via three dimensional Rossler oscillators and the macaque monkey data that was described in chapter 2. The results indicated a better performance of our Gaussian based PLV than that of the non- parametricapproach. The problem of cross-talk due to the poor spatial resolution is the limiting factor for the functional connectivity analysis in EEG and MEG data. The standard technique coherence usually leads to spurious connectivity results. Several approaches including IC, PLI, LC and OC have been proposed to overcome this problem. However, none of these techniques considers the cross-talk from the other sources in a multivariate network. Inchapter4,weprovideacomprehensivereviewoftheseapproachesandshow theinvariantpropertyofLCinpresenceofasinglecross-talk. Wethenproposeanovel measurePLCasanextensiontoLC.TheideabehindPLCistoexcludetheinstantaneous effects from the remaining signals using real coefficients and L1 penalty and compute the time-lagged interaction between the residuals. We evaulated the proposed approach as well as the other approaches using a realistic simulations of MEG data. The results indicatetherobustnessofPLCagainstmultiplecross-talksinamultivariatenetwork. Investigating functional connectivity networks from electrophysiological signals is still challenging despite the advanced technologies in neuroimaging. We need robust measures for a reliable estimation and interpretation. This dissertation is an attempt 98 to provide a better understanding to the existing measures and solutions to some well- knownproblemsinfunctionalconnectivitystudiesofelectrophysiologicaldata. 5.2 FutureDirections The findings in this dissertation can lead to some possible directions which requires further investigation. Below we describe possible future directions of the material pre- sentedinthisdissertation. 5.2.1 PartialPhaseSynchronization In chapter 3, we skipped the multivariate phase sychronization model proposed by Cadieu et al [CK10]. This statistical model extends the von Mises model introduced insection2.2.2tothemultivariatecasetocapturedirectphasecoupling. Theyintroduce anM-dimensionalphasedistributionforphasevectors p(jK) = 1 Z(K) exp [ 1 2 M ∑ m;n=1 mn cos(ϕ m ϕ n mn ) ] (5.1) = 1 Z(K) exp [ 1 2 x H Kx ] (5.2) where K mn = mn e im;n are the parameters of the distribution for all m and n, x m △ =e jϕm andZ(K)isthenormalizationconstantsothat ∫ 1 0 p(jK)d = 1. The parameters mn represent the strength of the direct coupling between nodes m and n and are nonzero only if m and n are coupled directly. The matrixK plays the same role as the concentration matrix (inverse covariance) in a multivariate Gaussian distribution. Thus, computation of this matrix allows us to estimate the direct phase coupling. Score matching method [Hyv06] is used to estimate parameters in equation 5.2whichhasnontrivialnormalizationconstant. 99 The partial PLV methods discussed in chapter 3 use second order statistics to sum- marizeinteractions. Cadieu’smethodratherdefinesamultivariateprobabilisticmodelto quantifyinteractions. TherelationshipbetweenCadieu’smethodandtheothermethods describedinchapter3isyettobeexploredindetail. 5.2.2 ConnectivityinthePresenceofCross-talk Inchapter4,weappliedconnectivitymeasurestoasimulateddataasweneededtoknow the gound truth to demonstrate the performances. However, all these measures have been proposed to be applicable to the real settings. Therefore, it would be interesting to see how these measures behave when applied to an experimental data. Definition of the ground truth would be challenging in this case. One approach would be to use the connectivity results from an intracranial data (e.g. ECoG) and compare the results with anoninvasivedataresultsfromthesamesubjectandtask(e.g. EEG,MEG).Inthiscase, definition of the region of interests would be crucial. They should be selected based on the locations of the intracranial electrodes. We will expect the PLC results to yield the similarnetworkresultsbetweenintracranialandnoninvasivedataifthereisacross-talk betweentheregions. 5.2.3 ClinicalApplications Anabnormalityinbrainfunctionalnetworkcanimplyphysiologicalandpsychiatricdis- eases such as dementia, schizophrenia, autism, attention deficit/hyperactivity disorder, depression and epilepsy [Men11]. Therefore, clinical application of the functional con- nectivity analysis may provide more insight for diagnosis and treatment. Furthermore, understanding functional connectivity in abnormal brains will improve the understand- ing of the dynamics of a healthy brain [HWA + 13]. These are some of the fundamental motivations behind the functional connectivity studies. We believe that the findings in 100 this dissertation will contribute to the functional connectivity studies of such diseases describedinsection1.4. 101 ReferenceList [ABES10] S. Aviyente, E.M. Bernat, W.S. Evans, and S.R. Sponheim. A phase synchrony measure for quantifying dynamic functional integration in thebrain. Human brainmapping,32(1):80–93,2010. [AC10] John JB Allen and Michael X Cohen. Deconstructing the “resting” state: exploring the temporal dynamics of frontal alpha asymmetry as anendophenotypefordepression. Frontiersinhumanneuroscience,4, 2010. [Ade87] GeorgeAdelman. Encyclopediaofneuroscience. 1987. [APL13] Sergul Aydore, Dimitrios Pantazis, and Richard M Leahy. A note on thephaselockingvalueanditsproperties. NeuroImage,2013. [ARN + 05] F. Amor, D. Rudrauf, V. Navarro, K. N’diaye, L. Garnero, J. Mar- tinerie, and M. Le Van Quyen. Imaging brain synchrony at high spatio-temporal resolution: application to meg signals during absence seizures. Signal processing,85(11):2101–2111,2005. [Bai12] Sylvain Baillet. Magnetoencephalography. In Brain Mapping, pages 77–89.Springer,2012. [BAK12] Gy¨ orgyBuzs´ aki,CostasAAnastassiou,andChristofKoch. Theorigin of extracellular fields and currents—eeg, ecog, lfp and spikes. Nature ReviewsNeuroscience,13(6):407–420,2012. [BCN93] S.L. Bressler, R. Coppola, and R. Nakamura. Episodic multiregional cortical coherence at multiple frequencies during visual task perfor- mance. Nature,366(6451):153–156,1993. [BDL + 04] A. Brovelli, M. Ding, A. Ledberg, Y. Chen, R. Nakamura, and S.L. Bressler. Beta oscillations in a large-scale sensorimotor cortical net- work: directional influences revealed by granger causality. Proceed- ingsoftheNationalAcademyofSciencesoftheUnitedStatesofAmer- ica,101(26):9849,2004. 102 [BDY99] S.L. Bressler, M. Ding, and W. Yang. Investigation of cooperative cortical dynamics by multivariate autoregressive modeling of event- relatedlocalfieldpotentials. Neurocomputing,26:625–631,1999. [BGdN06] Onureena Banerjee, Laurent El Ghaoui, Alexandre d’Aspremont, and Georges Natsoulis. Convex optimization techniques for fitting sparse gaussian graphical models. In Proceedings of the 23rd international conferenceon Machinelearning,pages89–96.ACM,2006. [BK68] HW Beams and RG Kessel. The golgi apparatus: structure and func- tion. Internationalreviewof cytology,23:209–276,1968. [BML01] S. Baillet, J.C. Mosher, and R.M. Leahy. Electromagnetic brain map- ping. Signal ProcessingMagazine,IEEE,18(6):14–30,2001. [BMLA + 06] Danielle S Bassett, Andreas Meyer-Lindenberg, Sophie Achard, Thomas Duke, and Edward Bullmore. Adaptive reconfiguration of fractal small-world human brain functional networks. Proceedings of the National Academy of Sciences,103(51):19518–19523,2006. [BNM + 12] Danielle S Bassett, Brent G Nelson, Bryon A Mueller, Jazmin Cam- chong, and Kelvin O Lim. Altered resting state complexity in schizophrenia. Neuroimage,59(3):2196–2207,2012. [BPFR01] J. Bhattacharya, H. Petsche, U. Feldmann, and B. Rescher. Eeg gamma-band phase synchronization between posterior and frontal cortex during mental rotation in humans. Neuroscience Letters, 311(1):29–32,2001. [BPS + 99] SN Baker, N Philbin, R Spinks, EM Pinches, DM Wolpert, DG Mac- Manus, Q Pauluis, and RN Lemon. Multiple single unit recording in thecortexofmonkeysusingindependentlymoveablemicroelectrodes. Journalof neurosciencemethods,94(1):5–17,1999. [Bre95] S.L. Bressler. Large-scale cortical networks and cognition. Brain ResearchReviews,20(3):288–304,1995. [BS09] Ed Bullmore and Olaf Sporns. Complex brain networks: graph theo- retical analysis of structural and functional systems. Nature Reviews Neuroscience,10(3):186–198,2009. [BST + 09] Randy L Buckner, Jorge Sepulcre, Tanveer Talukdar, Fenna M Krienen, Hesheng Liu, Trey Hedden, Jessica R Andrews-Hanna, Reisa A Sperling, and Keith A Johnson. Cortical hubs revealed by intrinsic functional connectivity: mapping, assessment of stability, 103 and relation to alzheimer’s disease. The Journal of Neuroscience, 29(6):1860–1873,2009. [BTVM05] Dale L Bailey, David W Townsend, Peter E Valk, and Michael N Maisey. Positronemission tomography. Springer,2005. [BWB + 01] Fabrice Bartolomei, Fabrice Wendling, Jean-Jacques Bellanger, Jean R´ egis, and Patrick Chauvel. Neural networks involving the medial temporal structures in temporal lobe epilepsy. Clinical neurophysiol- ogy,112(9):1746–1760,2001. [BWR + 04] F Bartolomei, F Wendling, J Regis, M Gavaret, M Guye, and P Chau- vel. Pre-ictalsynchronicityinlimbicnetworksofmesialtemporallobe epilepsy. Epilepsy research,61(1):89–104,2004. [BWV + 99] Fabrice Bartolomei, Fabrice Wendling, Jean-Pierre Vignal, Sylvia Kochen, Jean-Jacques Bellanger, Jean-Michel Badier, R´ egine Le Bouquin-Jeannes, and Patrick Chauvel. Seizures of temporal lobe epilepsy: identificationofsubtypesbycoherenceanalysisusingstereo- electro-encephalography. Clinical neurophysiology, 110(10):1741– 1754,1999. [BZYHH95] Bharat Biswal, F Zerrin Yetkin, Victor M Haughton, and James S Hyde. Functional connectivity in the motor cortex of resting human brain using echo-planar mri. Magnetic resonance in medicine, 34(4):537–541,1995. [cA] Course code AM 470701. Principles of neuroscience. Biomolecular Sciences MSc. [Car87] G Clifford Carter. Coherence and time delay estimation. Proceedings of the IEEE,75(2):236–255,1987. [CCK + 12] R.T. Canolty, C.F. Cadieu, K. Koepsell, K. Ganguly, R.T. Knight, and J.M. Carmena. Detecting event-related changes of multivariate phase coupling in dynamic brain networks. Journal of neurophysiology, 107(7):2020–2031,2012. [Cel07] P. Celka. Statistical analysis of the phase-locking value. Signal Pro- cessing Letters,IEEE,14(9):577–580,2007. [CFM + 10] Bruce Crosson, Anastasia Ford, Keith M McGregor, Marcus Meinzer, SergeyCheshkov,XiufengLi,DelainaWalker-Batson,andRichardW Briggs. Functional imaging and related techniques: An introduction for rehabilitation researchers. Journal of rehabilitation research and development,47(2):vii,2010. 104 [CHA + 00] Dietmar Cordes, Victor M Haughton, Konstantinos Arfanakis, Gary J Wendt, Patrick A Turski, Chad H Moritz, Michelle A Quigley, and MElizabethMeyerand. Mappingfunctionallyrelatedregionsofbrain withfunctionalconnectivitymrimaging. AmericanJournalofNeuro- radiology,21(9):1636–1644,2000. [CHA + 01] DietmarCordes,VictorMHaughton,KonstantinosArfanakis,JohnD Carew, Patrick A Turski, Chad H Moritz, Michelle A Quigley, and M Elizabeth Meyerand. Frequencies contributing to functional con- nectivityinthecerebralcortexin“resting-state”data. AmericanJour- nal of Neuroradiology,22(7):1326–1333,2001. [CJD + 03] Marco Catani, Derek K Jones, Rosario Donato, et al. Occipito- temporal connections in the human brain. Brain, 126(9):2093–2107, 2003. [CK91] RE Challis and RI Kitney. Biomedical signal processing (in four parts). Medical and Biological Engineering and Computing, 29(1):1– 17,1991. [CK10] C.F. Cadieu and K. Koepsell. Phase coupling estimation from multi- variatephasestatistics. Neuralcomputation,22(12):3107–3126,2010. [CLC + 99] Thomas E Conturo, Nicolas F Lori, Thomas S Cull, Erbil Akbudak, Abraham Z Snyder, Joshua S Shimony, Robert C McKinstry, Harold Burton,andMarcusERaichle.Trackingneuronalfiberpathwaysinthe livinghumanbrain. ProceedingsoftheNationalAcademyofSciences, 96(18):10422–10427,1999. [Dah00] R.Dahlhaus. Graphicalinteractionmodelsformultivariatetimeseries 1. Metrika,51(2):157–172,2000. [DBSB + 05] Robert F Dougherty, Michal Ben-Shachar, Roland Bammer, Alyssa A Brewer, and Brian A Wandell. Functional organization of human occipital-callosalfibertracts. ProceedingsoftheNationalAcademyof Sciences of the United States of America,102(20):7350–7355,2005. [DBYL00] M.Ding,S.L.Bressler,W.Yang,andH.Liang. Short-windowspectral analysis of cortical event-related potentials by adaptive multivariate autoregressive modeling: data preprocessing, model validation, and variabilityassessment. Biologicalcybernetics,83(1):35–45,2000. [DCB91] Henri M Duvernoy, Emmanuel A Cabanis, and P Bourgouin. The humanbrain: surface,three-dimensionalsectionalanatomyandMRI. SpringerVerlagWien,1991. 105 [DCS99] Alain Destexhe, Diego Contreras, and Mircea Steriade. Spatiotempo- ral analysis of local field potentials and unit discharges in cat cerebral cortex during natural wake and sleep states. The Journal of Neuro- science,19(11):4595–4608,1999. [DGCV02] O. David, L. Garnero, D. Cosmelli, and F.J. Varela. Estimation of neuraldynamicsfrommeg/eegcorticalcurrentdensitymaps: applica- tiontothereconstructionoflarge-scalecorticalsynchrony. Biomedical Engineering,IEEE Transactionson,49(9):975–987,2002. [DLF + 00] Anders M Dale, Arthur K Liu, Bruce R Fischl, Randy L Buckner, John W Belliveau, Jeffrey D Lewine, and Eric Halgren. Dynamic statistical parametric mapping: combining fmri and meg for high- resolutionimagingofcorticalactivity. Neuron,26(1):55–67,2000. [DLSDS + 05] Marilena De Luca, Stephen Smith, Nicola De Stefano, Antonio Fed- erico, and Paul M Matthews. Blood oxygenation level dependent contrast resting state networks are relevant to functional activity in the neocortical sensorimotor system. Experimental Brain Research, 167(4):587–594,2005. [DR58] W.B. Davenport and W.L. Root. An introduction to the theory of ran- dom signals and noise,volume11. McGraw-HillNewYork,1958. [DRKW08] S.M. Doesburg, A.B. Roggeveen, K. Kitajo, and L.M. Ward. Large- scale gamma-band phase synchronization and selective attention. CerebralCortex,18(2):386–396,2008. [DS92] Robert B Duckrow and Susan S Spencer. Regional coherence and the transfer of ictal activity during seizure onset in the medial tem- poral lobe. Electroencephalography and clinical neurophysiology, 82(6):415–422,1992. [DS93] Anders M Dale and Martin I Sereno. Improved localizadon of cor- tical activity by combining eeg and meg with mri cortical surface reconstruction: A linear approach. Journal of cognitive neuroscience, 5(2):162–176,1993. [EMFO05] AndreasKEngel,ChristianKEMoll,ItzhakFried,andGeorgeAOje- mann. Invasiverecordingsfromthehumanbrain: clinicalinsightsand beyond. NatureReviewsNeuroscience,6(1):35–47,2005. [EMZ + 12] Arne Ewald, Laura Marzetti, Filippo Zappasodi, Frank C Meinecke, and Guido Nolte. Estimating true brain connectivity from eeg/meg 106 datainvarianttolinearandstatictransformationsinsensorspace. Neu- roImage,60(1):476–488,2012. [ETS + 13] JeromeEngel,PaulMThompson,JohnMStern,RichardJStaba,Ana- tol Bragin, Istvan Mody, et al. Connectomics and epilepsy. Current opinion in neurology,26(2):186–194,2013. [FCS + 06] Michael D Fox, Maurizio Corbetta, Abraham Z Snyder, Justin L Vin- cent, and Marcus E Raichle. Spontaneous neuronal activity distin- guishes human dorsal and ventral attention systems. Proceedings of the National Academy of Sciences,103(26):10046–10051,2006. [FET00] Thomas C Ferree, K Jeffrey Eriksen, and Don M Tucker. Regional head tissue conductivity estimation for improved eeg analysis. Biomedical Engineering, IEEE Transactions on, 47(12):1584–1592, 2000. [FHT08] JeromeFriedman,TrevorHastie,andRobertTibshirani.Sparseinverse covarianceestimationwiththegraphicallasso. Biostatistics,9(3):432– 441,2008. [FR07] Michael D Fox and Marcus E Raichle. Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. NatureReviewsNeuroscience,8(9):700–711,2007. [FRBM88] George Fein, Jonathan Raz, Fiona F Brown, and Edward L Mer- rin. Common reference coherence data are confounded by power and phase effects. Electroencephalography and clinical neurophysiology, 69(6):581–584,1988. [Fri09] Karl Friston. Causal modelling and brain connectivity in functional magneticresonanceimaging. PLoS biology,7(2):e1000033,2009. [Fri11] Karl J Friston. Functional and effective connectivity: a review. Brain connectivity,1(1):13–36,2011. [FSZR05] Michael D Fox, Abraham Z Snyder, Jeffrey M Zacks, and Mar- cusERaichle. Coherentspontaneousactivityaccountsfortrial-to-trial variability in human evoked brain responses. Nature neuroscience, 9(1):23–25,2005. [GA06] Yosef Grodzinsky and Katrin Amunts. Broca’s region. Oxford Uni- versityPressNewYork,2006. [Gal08] R.G. Gallager. Circularly-symmetric gaussian random vectors. preprint,pages1–9,2008. 107 [GBR08] MaximeGuye,FabriceBartolomei,andJean-PhilippeRanjeva. Imag- ing structural and functional connectivity: towards a unified defini- tion of human brain organization? Current opinion in neurology, 21(4):393–403,2008. [GDMH + 00] S Goncalves, JC De Munck, RM Heethaar, FH Lopes Da Silva, and BWVanDijk. Theapplicationofelectricalimpedancetomographyto reducesystematicerrorsintheeeginverseproblem-asimulationstudy. Physiologicalmeasurement,21(3):379,2000. [GG01] Murray Grossman and Jean Gotman. As time goes by high temporal and spatial resolution in cognitively related cortical function. Neurol- ogy,57(11):1947–1948,2001. [GHAEC08] Germ´ an G´ omez-Herrero, Mercedes Atienza, Karen Egiazarian, and Jose L Cantero. Measuring directional coupling between eeg sources. Neuroimage,43(3):497–508,2008. [GL96] J Gotman and V Levtova. Amygdala-hippocampus relationships in temporal lobe seizures: a phase-coherence study. Epilepsy research, 25(1):51–57,1996. [Glo85] Pierre Gloor. Neuronal generators and the problem of localiza- tion in electroencephalography: application of volume conductor the- ory to electroencephalography. Journal of clinical neurophysiology, 2(4):327–354,1985. [GMDV + 11] Juli´ anJGonz´ alez,SoledadMa˜ nas,Lu´ ısDeVera,LeopoldoDM´ endez, Santiago L´ opez, Jos´ e M Garrido, and Ernesto Pereda. Assessment of electroencephalographic functional connectivity in term and preterm neonates. ClinicalNeurophysiology,122(4):696–702,2011. [Goo63] NRGoodman. Statisticalanalysisbasedonacertainmultivariatecom- plex gaussian distribution (an introduction). The Annals of mathemat- ical statistics,34(1):152–177,1963. [GPO + 10] Alexandre Gramfort, Th´ eodore Papadopoulo, Emmanuel Olivi, Mau- reen Clerc, et al. Openmeeg: opensource software for quasistatic bio- electromagnetics. Biomed. Eng.Online,9(1):45,2010. [Gra07] Gunnar Grant. How the 1906 nobel prize in physiology or medicine was shared between golgi and cajal. Brain research reviews, 55(2):490–498,2007. 108 [Gre08] Michael Greicius. Resting-state functional connectivity in neuropsy- chiatricdisorders.Currentopinioninneurology,21(4):424–430,2008. [GVN + 05] R. Guevara, J.L.P. Velazquez, V. Nenadovic, R. Wennberg, G. Sen- janovi´ c, and L.G. Dominguez. Phase synchronization measure- ments using electroencephalographic recordings. Neuroinformatics, 3(4):301–313,2005. [HBB + 12] Arjan Hillebrand, Gareth R Barnes, Johannes L Bosboom, Henk W Berendse, and Cornelis J Stam. Frequency-dependent functional con- nectivity within resting-state networks: an atlas-based meg beam- formersolution. Neuroimage,59(4):3909–3921,2012. [HCE07] Yong He, Zhang J Chen, and Alan C Evans. Small-world anatomical networks in the human brain revealed by cortical thickness from mri. Cerebralcortex,17(10):2407–2419,2007. [HCT + 10] JPaulHamilton,GangChen,MoriahEThomason,MirraESchwartz, and Ian H Gotlib. Investigating neural primacy in major depressive disorder: multivariate granger causality analysis of resting-state fmri time-seriesdata. Molecularpsychiatry,16(7):763–772,2010. [HE10] Yong He and Alan Evans. Graph theoretical modeling of brain con- nectivity. Currentopinion in neurology,23(4):341–350,2010. [HHC + 12] Joerg F Hipp, David J Hawellek, Maurizio Corbetta, Markus Siegel, and Andreas K Engel. Large-scale cortical correlation structure of spontaneousoscillatoryactivity. Natureneuroscience,15(6):884–890, 2012. [HHI + 93] MattiH¨ am¨ al¨ ainen,RiittaHari,RistoJIlmoniemi,JukkaKnuutila,and Olli V Lounasmaa. Magnetoencephalography—theory, instrumenta- tion, and applications to noninvasive studies of the working human brain. Reviewsof modern Physics,65(2):413,1993. [HL08] Bin He and Zhongming Liu. Multimodal functional neuroimaging: integrating functional mri and eeg/meg. Biomedical Engineering, IEEE Reviewsin,1:23–40,2008. [HM77] DanielBHierandJPMohr.Incongruousoralandwrittennaming: Evi- dence for a subdivision of the syndrome of wernicke’s aphasia. Brain and language,4(1):115–126,1977. [HML99] MX Huang, JC Mosher, and RM Leahy. A sensor-weighted overlapping-sphere head model and exhaustive head model compar- isonformeg. Physicsin medicine and biology,44(2):423,1999. 109 [HMO + 87] Bin He, TOSHIMITSU Musha, YOSHIWO Okamoto, Saburo Homma, Yoshio Nakajima, and Toshio Sato. Electric dipole tracing in the brain by means of the boundary element method and its accu- racy. Biomedical Engineering, IEEE Transactions on, (6):406–414, 1987. [Hor03] B. Horwitz. The elusive concept of brain connectivity. Neuroimage, 19(2):466–470,2003. [HPBL10] Hua Brian Hui, Dimitrios Pantazis, Steven L Bressler, and Richard M Leahy. Identifying true cortical interactions in meg using the nulling beamformer. NeuroImage,49(4):3161–3174,2010. [HRPH10] AapoHyv¨ arinen,PavanRamkumar,LauriParkkonen,andRiittaHari. Independent component analysis of short-time fourier transforms for spontaneouseeg/meganalysis. Neuroimage,49(1):257–271,2010. [HRS04] J.M.Hurtado,L.L.Rubchinsky,andK.A.Sigvardt. Statisticalmethod fordetectionofphase-lockingepisodesinneuraloscillations. Journal of neurophysiology,91(4):1883–1898,2004. [HS89] MATTISHamalainenandJukkaSarvas. Realisticconductivitygeom- etrymodelofthehumanheadforinterpretationofneuromagneticdata. BiomedicalEngineering,IEEETransactionson,36(2):165–171,1989. [HSM04] Scott A Huettel, Allen W Song, and Gregory McCarthy. Functional magnetic resonance imaging, volume 1. Sinauer Associates Sunder- land,MA,2004. [HSZ + 08] Biyu J He, Abraham Z Snyder, John M Zempel, Matthew D Smyth, and Marcus E Raichle. Electrophysiological correlates of the brain’s intrinsic large-scale functional architecture. Proceedings of the National Academy of Sciences,105(41):16039–16044,2008. [HTF + 09] TrevorHastie,RobertTibshirani,JeromeFriedman, THastie,JFried- man,andRTibshirani. Theelementsofstatisticallearning,volume2. Springer,2009. [HTM99] BarryHorwitz,MATagamets,andAnthonyRandalMcIntosh. Neural modeling,functionalbrainimaging,andcognition. Trendsincognitive sciences,3(3):91–98,1999. [HWA + 13] R Matthew Hutchison, Thilo Womelsdorf, Elena A Allen, Peter A Bandettini, Vince D Calhoun, Maurizio Corbetta, Stefania 110 Della Penna, Jeff H Duyn, Gary H Glover, Javier Gonzalez- Castillo, et al. Dynamic functional connectivity: promise, issues, and interpretations. Neuroimage,80:360–378,2013. [HYS + 10] DorotheaIHorn,ChunshuiYu,JohannSteiner,JuliaBuchmann,Joern Kaufmann, Annemarie Osoba, Ulf Eckert, Kathrin C Zierhut, Kolja Schiltz, Huiguang He, et al. Glutamatergic and resting-state func- tional connectivity correlates of severity in major depression–the role ofpregenualanteriorcingulatecortexandanteriorinsula. Frontiersin systems neuroscience,4,2010. [Hyv06] A. Hyvarinen. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(1):695, 2006. [IMCRMG + 07] Y Iturria-Medina, EJ Canales-Rodriguez, L Melie-Garcia, PA Valdes- Hernandez, E Martinez-Montes, Y Aleman-Gomez, and JM S´ anchez- Bornot. Characterizing brain anatomical connections using diffusion weightedmriandgraphtheory. Neuroimage,36(3):645–660,2007. [JP70] EGJonesandTPSPowell. Ananatomicalstudyofconvergingsensory pathwayswithinthecerebralcortexofthemonkey. Brain,93(4):793– 820,1970. [JS01] S.R. Jammalamadaka and A. Sengupta. Topics in circular statistics, volume5. WorldScientificPubCoInc,2001. [JVM + 12] David T Jones, Prashanthi Vemuri, Matthew C Murphy, Jeffrey L Gunter, Matthew L Senjem, Mary M Machulda, Scott A Przybel- ski, Brian E Gregg, Kejal Kantarci, David S Knopman, et al. Non- stationarity in the “resting brain’s” modular architecture. PloS one, 7(6):e39731,2012. [KDTB01] Maciej Kami´ nski, Mingzhou Ding, Wilson A Truccolo, and Steven L Bressler. Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of sig- nificance. Biologicalcybernetics,85(2):145–157,2001. [KKB04] Rafal Kus, Maciej Kaminski, and Katarzyna J Blinowska. Determi- nationofeegactivitypropagation: pair-wiseversusmultichannelesti- mate. Biomedical Engineering, IEEE Transactions on, 51(9):1501– 1510,2004. 111 [KNB + 09] Steffen Katzner, Ian Nauhaus, Andrea Benucci, Vincent Bonin, Dario L Ringach, and Matteo Carandini. Local origin of field poten- tialsinvisualcortex. Neuron,61(1):35–41,2009. [Koc08] H Koch. Squids–quo vadis? Physica C: Superconductivity, 468(15):1112–1114,2008. [KSJS06] A. Klein, T. Sauer, A. Jedynak, and W. Skrandies. Conventional and wavelet coherence applied to sensory-evoked electrical brain activity. BiomedicalEngineering,IEEETransactionson,53(2):266–272,2006. [Lau96] Steffen L Lauritzen. Graphical models. Oxford University Press, 1996. [LB03] Denis Le Bihan. Looking into the functional architecture of the brain withdiffusionmri.NatureReviewsNeuroscience,4(6):469–480,2003. [LCL03] J.P.Lachaux,M.Chavez,andA.Lutz.Asimplemeasureofcorrelation across time, frequency and space between continuous brain signals. Journalof neurosciencemethods,123(2):175–188,2003. [LDB01] H. Liang, M. Ding, and S.L. Bressler. Temporal dynamics of infor- mation flow in the cerebral cortex. Neurocomputing, 38:1429–1435, 2001. [LFdZD10] Zhongming Liu, Masaki Fukunaga, Jacco A de Zwart, and Jeff H Duyn. Large-scale spontaneous fluctuations and correlations in brain electrical activity observed with magnetoencephalography. Neuroim- age,51(1):102–111,2010. [LH08] Zhongming Liu and Bin He. fmri–eeg integrated cortical source imaging by use of time-variant spatial constraints. NeuroImage, 39(3):1198–1214,2008. [LLG06] Frederick JP Langheim, Arthur C Leuthold, and Apostolos P Geor- gopoulos. Synchronous dynamic brain networks revealed by magne- toencephalography. Proceedings of the National Academy of Sciences of the United States of America,103(2):455–459,2006. [LMS98] MJLowe,BJMock,andJASorenson. Functionalconnectivityinsin- gleandmultisliceechoplanarimagingusingresting-statefluctuations. Neuroimage,7(2):119–132,1998. [LPA + 01] Nikos K Logothetis, Jon Pauls, Mark Augath, Torsten Trinath, and AxelOeltermann. Neurophysiologicalinvestigationofthebasisofthe fmrisignal. Nature,412(6843):150–157,2001. 112 [LRM + 99] J.P.Lachaux,E.Rodriguez,J.Martinerie,F.J.Varela,etal. Measuring phase synchrony in brain signals. Human brain mapping, 8(4):194– 208,1999. [LRM + 00] J.P. Lachaux, E. Rodriguez, L.E.V.A.N.Q. MICHEL, L. ANTOINE, J. Martinerie, and J.V. FRANCISCO. Studying single-trials of phase synchronousactivityinthebrain. InternationalJournalofBifurcation and Chaos,10(10):2429–2439,2000. [LS05] H.O.LancasterandE.Seneta. Chi-SquareDistribution. WileyOnline Library,2005. [Luc17] Keith Lucas. The conduction of the nervous impulse. Longmans, GreenandCompany,1917. [LVQFL + 01] M. Le Van Quyen, J. Foucher, J.P. Lachaux, E. Rodriguez, A. Lutz, J. Martinerie, and F.J. Varela. Comparison of hilbert transform and wavelet methods for the analysis of neuronal synchrony. Journal of neurosciencemethods,111(2):83–98,2001. [Mac02] Malcolm Macmillan. An odd kind of fame: Stories of Phineas Gage. MITPress,2002. [MBM + 10] Francesco Musso, J¨ urgen Brinkmeyer, Arian Mobascher, Tracy War- brick, and Georg Winterer. Spontaneous brain activity and eeg microstates.anoveleeg/fmrianalysisapproachtoexploreresting-state networks. Neuroimage,52(4):1149–1161,2010. [MDGN08] Laura Marzetti, Cosimo Del Gratta, and Guido Nolte. Understanding brain connectivity from eeg data by identifying systems composed of interactingsources. Neuroimage,42(1):87–98,2008. [Men11] Vinod Menon. Large-scale brain networks and psychopathology: a unifying triple network model. Trends in cognitive sciences, 15(10):483–506,2011. [MLDEE00] F.Mormann,K.Lehnertz,P.David,andC.EElger. Meanphasecoher- ence as a measure for phase synchronization and its application to the eeg of epilepsy patients. Physica D: Nonlinear Phenomena, 144(3- 4):358–369,2000. [MLL92] John C Mosher, Paul S Lewis, and Richard M Leahy. Multiple dipole modelingandlocalizationfromspatio-temporalmegdata. Biomedical Engineering,IEEE Transactionson,39(6):541–557,1992. 113 [MM04] C.J. Mecklin and D.J. Mundfrom. An appraisal and bibliography of tests for multivariate normality. International Statistical Review, 72(1):123–138,2004. [MMH00] T. Mima, T. Matsuoka, and M. Hallett. Functional coupling of human rightandleftcorticalmotorareasdemonstratedwithpartialcoherence analysis. Neuroscienceletters,287(2):93–96,2000. [MP95] Jaakko Malmivuo and Robert Plonsey. Bioelectromagnetism: princi- ples and applications of bioelectric and biomagnetic fields. Oxford UniversityPress,1995. [MPDG + 07] DanteMantini,MauroGianniPerrucci,CosimoDelGratta,GianLuca Romani, and Maurizio Corbetta. Electrophysiological signatures of restingstatenetworksinthehumanbrain. ProceedingsoftheNational Academy of Sciences,104(32):13170–13175,2007. [NBT + 05] Katherine L Narr, Robert M Bilder, Arthur W Toga, Roger P Woods, DavidERex,PhilipRSzeszko,DelbertRobinson,SergeSevy,Handan Gunduz-Bruce, Yung-Ping Wang, et al. Mapping cortical thickness andgraymatterconcentrationinfirstepisodeschizophrenia. Cerebral Cortex,15(6):708–719,2005. [NBW + 04] G. Nolte, O. Bai, L. Wheaton, Z. Mari, S. Vorbach, M. Hallett, et al. Identifying true brain interaction from eeg data using the imaginary part of coherency. Clinical Neurophysiology, 115(10):2292–2307, 2004. [NMD + 08] Yuval Nir, Roy Mukamel, Ilan Dinstein, Eran Privman, Michal Harel, Lior Fisch, Hagar Gelbard-Sagiv, Svetlana Kipervasser, Fani Andel- man, Miri Y Neufeld, et al. Interhemispheric correlations of slow spontaneous neuronal fluctuations revealed in human sensory cortex. Natureneuroscience,11(9):1100–1108,2008. [NSS + 99] Paul L Nunez, Richard B Silberstein, Zhiping Shi, Matthew R Car- penter, Ramesh Srinivasan, Don M Tucker, Scott M Doran, Peter J Cadusch, and Ranjith S Wijesinghe. Eeg coherency ii: experimen- tal comparisons of multiple measures. Clinical Neurophysiology, 110(3):469–486,1999. [NSW + 97] P.L. Nunez, R. Srinivasan, A.F. Westdorp, R.S. Wijesinghe, D.M. Tucker, R.B. Silberstein, and P.J. Cadusch. Eeg coherency:: I: statis- tics,referenceelectrode,volumeconduction,laplacians,corticalimag- ing,andinterpretationatmultiplescales. Electroencephalographyand clinical Neurophysiology,103(5):499–515,1997. 114 [Nun95] Paul L Nunez. Neocortical dynamics and human EEG rhythms. OxfordUniversityPress,USA,1995. [Nun06] Paul L Nunez. Electric fields of the brain: the neurophysics of EEG. OxfordUniversityPress,2006. [NZN + 08] Guido Nolte, Andreas Ziehe, Vadim V Nikulin, Alois Schl¨ ogl, Nicole Kr¨ amer, Tom Brismar, and Klaus-Robert M¨ uller. Robustly estimating theflowdirectionofinformationincomplexphysicalsystems. Physi- cal reviewletters,100(23):234101,2008. [OGT + 10] A. Ossadtchi, RE Greenblatt, VL Towle, MH Kohrman, and K. Kamada. Inferring spatiotemporal network patterns from intracra- nialeegdata. ClinicalNeurophysiology,121(6):823–835,2010. [PC36] Linus Pauling and Charles D Coryell. The magnetic properties and structure of hemoglobin, oxyhemoglobin and carbonmonoxyhe- moglobin. Proceedings of the National Academy of Sciences of the United States of America,22(4):210,1936. [Pev02] Jonathan Pevsner. Leonardo da vinci’s contributions to neuroscience. TRENDS in Neurosciences,25(4):217–220,2002. [PGH + 03] Gert Pfurtscheller, Bernard Graimann, Jane E Huggins, Simon P Levine, and Lori A Schuh. Spatiotemporal patterns of beta desyn- chronizationandgammasynchronizationincorticographicdataduring self-paced movement. Clinical Neurophysiology, 114(7):1226–1236, 2003. [Plo69] RobertPlonsey. Bioelectric phenomena. WileyOnlineLibrary,1969. [PM07] Roberto D Pascual-Marqui. Coherence and phase synchronization: generalizationtopairsofmultivariatetimeseries,andremovalofzero- lagcontributions. arXiv preprintarXiv:0706.1776,2007. [QKKG02] R.Q.Quiroga,A.Kraskov,T.Kreuz,andP.Grassberger. Performance of different synchronization measures in real data: a case study on electroencephalographic signals. Physical Review E, 65(4):041903, 2002. [RDK + 06] D. Rudrauf, A. Douiri, C. Kovach, J.P. Lachaux, D. Cosmelli, M.Chavez,C.Adam,B.Renault,J.Martinerie,andM.LeVanQuyen. Frequency flows and the time-frequency dynamics of multivariate phase synchronization in brain signals. Neuroimage, 31(1):209–227, 2006. 115 [RM06] MarcusERaichleandMarkAMintun. Brainworkandbrainimaging. Annu. Rev.Neurosci.,29:449–476,2006. [Rou99] Mark S Roulston. Estimating the errors on measured entropy and mutual information. Physica D: Nonlinear Phenomena, 125(3):285– 294,1999. [SCWM06] Andrew B Schwartz, X Tracy Cui, Douglas J Weber, and Daniel W Moran. Brain-controlledinterfaces: movementrestorationwithneural prosthetics. Neuron,52(1):205–220,2006. [SDMGIM + 10] Gretel Sanabria-Diaz, Lester Melie-Garc´ ıa, Yasser Iturria-Medina, Yasser Alem´ an-G´ omez, Gertrudis Hern´ andez-Gonz´ alez, Lourdes Vald´ es-Urrutia,L´ ıdiceGal´ an,andPedroVald´ es-Sosa.Surfaceareaand corticalthicknessdescriptorsrevealdifferentattributesofthestructural humanbrainnetworks. Neuroimage,50(4):1497–1510,2010. [SG09] Jan-MathijsSchoffelenandJoachimGross.Sourceconnectivityanaly- siswithmegandeeg. Humanbrainmapping,30(6):1857–1865,2009. [She91] Gordon M Shepherd. Foundations of the neuron doctrine. Oxford UnivPress,1991. [Sil95] RB Silberstein. Neocortical dynamics and human eeg rhythms. Neo- cortical Dynamics and Human EEG Rhythms,page708,1995. [SL08] Amir Shmuel and David A Leopold. Neuronal correlates of spon- taneous fluctuations in fmri signals in monkey visual cortex: impli- cations for functional connectivity at rest. Human brain mapping, 29(7):751–761,2008. [SMR + 08] Kaustubh Supekar, Vinod Menon, Daniel Rubin, Mark Musen, and Michael D Greicius. Network analysis of intrinsic functional brain connectivity in alzheimer’s disease. PLoS computational biology, 4(6):e1000100,2008. [SND07] C.J. Stam, G. Nolte, and A. Daffertshofer. Phase lag index: assess- ment of functional connectivity from multi channel eeg and meg with diminished bias from common sources. Human brain mapping, 28(11):1178–1193,2007. [SP98] JGStinstraandMJPeters. Thevolumeconductormayactasatempo- ral filter on the ecg and eeg. Medical and Biological Engineering and Computing,36(6):711–716,1998. 116 [SPK + 10] ¨ Unal Sako˘ glu, Godfrey D Pearlson, Kent A Kiehl, Y Michelle Wang, Andrew M Michael, and Vince D Calhoun. A method for evaluating dynamic functional network connectivity and task-modulation: appli- cation to schizophrenia. Magnetic Resonance Materials in Physics, Biologyand Medicine,23(5-6):351–366,2010. [SRM + 07] Christian Sorg, Valentin Riedl, Mark M¨ uhlau, Vince D Calhoun, Tom Eichele, Leonhard L¨ aer, Alexander Drzezga, Hans F¨ orstl, Alexander Kurz, Claus Zimmer, et al. Selective changes of resting-state net- works in individuals at risk for alzheimer’s disease. Proceedings of the National Academy of Sciences,104(47):18760–18765,2007. [STK05] OlafSporns,GiulioTononi,andRolfK¨ otter. Thehumanconnectome: astructural description ofthe human brain. PLoS computational biol- ogy,1(4):e42,2005. [SWD + 06] B. Schelter, M. Winterhalder, R. Dahlhaus, J. Kurths, and J. Timmer. Partial phase synchronization for multivariate synchronizing systems. Physical reviewletters,96(20):208103,2006. [Swe96] J.A. Swets. Signal detection theory and ROC analysis in psychology and diagnostics: Collected papers. Lawrence Erlbaum Associates, Inc,1996. [TBM + 11] Franc ¸ois Tadel, Sylvain Baillet, John C Mosher, Dimitrios Pantazis, and Richard M Leahy. Brainstorm: a user-friendly application for meg/eeg analysis. Computational intelligence and neuroscience, 2011:8,2011. [TC69] Raymond Carl Truex and Malcolm B Carpenter. Human neu- roanatomy. Williams&Wilkins,1969. [TE98] G. Tononi and G.M. Edelman. Consciousness and complexity. Sci- ence,282(5395):1846–1851,1998. [Tib96] Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. Series B (Methodological), pages267–288,1996. [TRW + 98] P. Tass, MG Rosenblum, J. Weule, J. Kurths, A. Pikovsky, J. Volk- mann,A.Schnitzler,andH.J.Freund. Detectionofn: mphaselocking from noisy data: application to magnetoencephalography. Physical ReviewLetters,81(15):3291–3294,1998. 117 [VL92] J Van Laere. [vesalius and the nervous system]. Verhandelingen- Koninklijke Academie voor Geneeskunde van Belgie, 55(6):533–576, 1992. [VLRM01] F.Varela,J.P.Lachaux,E.Rodriguez,andJ.Martinerie.Thebrainweb: phasesynchronizationandlarge-scaleintegration. Naturereviewsneu- roscience,2(4):229–239,2001. [VON11] A.Vo,S.Oraintara,andN.Nguyen.Vonndistributionofrelativephase forstatisticalimagemodelingincomplexwaveletdomain. SignalPro- cessing,91(1):114–125,2011. [VOvW + 11] M.Vinck,R.Oostenveld,M.vanWingerden,F.Battaglia,andC.Pen- nartz. Animprovedindexofphase-synchronizationforelectrophysio- logical data in the presence of volume-conduction, noise and sample- sizebias. Neuroimage,55(4):1548–1565,2011. [VVWW + 10] M. Vinck, M. Van Wingerden, T. Womelsdorf, P. Fries, and C. Pen- nartz.Thepairwisephaseconsistency: abias-freemeasureofrhythmic neuronalsynchronization. Neuroimage,51(1):112–122,2010. [Whi09] J. Whittaker. Graphical models in applied multivariate statistics. WileyPublishing,2009. [XPGF99] Jinhu Xiong, Lawrence M Parsons, Jia-Hong Gao, and Peter T Fox. Interregional connectivity to primary motor cortex revealed using mri restingstateimages. Humanbrainmapping,8(2-3):151–156,1999. 118 AppendixA PhaseSynchronization A.1 BiasofSquaredPLV Herewederiveanexpressionforthebiasof ^ PLV 2 sample . E [ ^ PLV 2 sample ] = E 2 4 1 N N ∑ n=1 e j∆ϕn 2 3 5 = 1 N 2 E [ N ∑ n=1 e j∆ϕn N ∑ m=1 e j∆ϕm ] (A.1) = 1 N 2 ( N + N ∑ n=1 N ∑ m̸=n E [ e j(∆ϕn∆ϕm) ] ) (A.2) Assuming ∆ϕ n and ∆ϕ m areindependentforn̸=m = 1 N 2 ( N + N ∑ n=1 N ∑ m̸=n E [ e j∆ϕn ] E [ e j∆ϕm ] ) (A.3) = 1 N 2 ( N +(N 2 N) E [ e j∆ϕ ] 2 ) = 1 N + ( 1 1 N ) PLV 2 (A.4) A.2 UnbiasedsquarePLVandPPC Rearrangingequ. (A.4),givesthefollowingunbiasedestimateofPLV 2 : ^ PLV 2 ub sample = 1 N1 ( ^ PLV 2 sample N1 ) (A.5) 119 Further manipulation of this expression shows it is equivalent to the pairwise phase consistency(PPC)measureofVinck. = 1 N1 2 4 1 N N ∑ n=1 e j∆ϕn 2 1 3 5 (A.6) = 1 N(N1) [ N + N ∑ n=1 N ∑ m̸=n e j(∆ϕm∆ϕn) ] 1 N1 (A.7) = 1 N(N1) N ∑ n=1 N ∑ m̸=n e j(∆ϕm∆ϕn) | {z } ∑ N1 n=1 ∑ N m>n e j(∆ϕm∆ϕn) + ∑ N n=1 ∑ N1 m<n e j(∆ϕm∆ϕn) (A.8) = 2 N(N1) N1 ∑ n=1 N ∑ m=n+1 cos(∆ϕ m ∆ϕ n ) (A.9) whichisthePPCmeasureinequ. 2.12in[VVWW + 10]. A.3 Conditional Distribution of Relative Phase for the BivariateCircularlySymmetricComplexGaussian First,weneedtofindanexpressionforthemarginaldistributionofamplitudesA(t) p(A(t)) = ∫ 2 0 ∫ 2 0 p(A(t);(t))dϕ 1 (t)dϕ 2 (t) (A.10) fromequ. (2.19) = A 1 (t)A 2 (t) 2 jK z j exp ( 2 ∑ i=1 ii A 2 i (t) ) ∫ 2 0 ∫ 2 0 exp(2 12 A 1 (t)A 2 (t)cos(ϕ 1 (t)ϕ 2 (t) 12 ))dϕ 1 (t)dϕ 2 (t) | {z } 4 2 I0(212A1(t)A2(t)) (A.11) = 4A 1 (t)A 2 (t)I 0 (2 12 A 1 (t)A 2 (t)) jK z j exp ( 2 ∑ i=1 ii A 2 i (t) ) : (A.12) 120 Substitutingthisresultintoequ. (2.20),weget p((t)jA(t)) = A 1 (t)A 2 (t) 2 jKzj exp { ∑ 2 i=1 ii A 2 i (t)2 12 A 1 (t)A 2 (t)cos(ϕ 1 (t)ϕ 2 (t) 12 ) } 4A 1 (t)A 2 (t)I 0 (2 12 A 1 (t)A 2 (t)) jKzj exp ( ∑ 2 i=1 ii A 2 i (t) ) = 1 4 2 I 0 (2 12 A 1 (t)A 2 (t)) e 2 12 A 1 (t)A 2 (t)cos(ϕ 1 (t)ϕ 2 (t) 12 ) (A.13) Then, the distribution of relative phase ∆ϕ(t) conditioned on the amplitudesA can be writtenas p(∆ϕ(t)jA(t)) = 1 2I 0 (2 12 A 1 (t)A 2 (t)) e 2 12 A 1 (t)A 2 (t)cos(∆ϕ(t) 12 ) : (A.14) A.4 The PLV for the Bivariate Circularly Symmetric ComplexGaussian E[e j∆ϕ(t) ] = E [ E [ e j∆ϕ(t) jA 1 (t);A 2 (t) ]] (A.15) usingthefirstmomentgeneratingfunctionofvonMisesdistribution: = E [ I 1 (2 12 A 1 (t)A 2 (t)) I 0 (2 12 A 1 (t)A 2 (t)) ] (A.16) = ∫ 1 0 ∫ 1 0 I 1 (2 12 A 1 (t)A 2 (t)) I 0 (2 12 A 1 (t)A 2 (t)) p(A 1 (t);A 2 (t))dA 1 (t)dA 2 (t) (A.17) 121 usingtheresultinequ. A.12: = 4 det(K z ) ∫ 1 0 ∫ 1 0 I 1 (2 12 A 1 (t)A 2 (t)) | {z } 1 ∫ =2 =2 e 2 12 A 1 (t)A 2 (t)sin() sin()d A 1 (t)A 2 (t)e 11 A 2 1 (t) 22 A 2 2 (t) dA 1 (t)dA 2 (t) (A.18) = 4 det(K z ) ∫ 1 0 ∫ 1 0 ∫ =2 =2 A 1 (t)A 2 (t)sin() e 11 A 2 1 (t)+2 12 A 1 (t)A 2 (t)sin 22 A 2 2 (t) dA 1 (t)dA 2 (t)d (A.19) = 4 det(K z ) ∫ =2 =2 sin() ∫ 1 0 A 2 (t)e 22 A 2 2 (t) ∫ 1 0 A 1 (t)e 11 A 2 1 (t)+2 12 A 1 (t)A 2 (t)sin() dA 1 (t) | {z } ≜B(;A 2 (t); 11 ; 12 ) dA 2 (t)d (A.20) B(;A 2 (t); 11 ; 12 ) = ∫ 1 0 A 1 (t)e 11 ( A 2 1 (t) 2 12 A 2 (t)sin() 11 A 1 (t) ) dA 1 (t) (A.21) = e 2 12 A 2 2 (t)sin 2 () 11 ∫ 1 0 A 1 (t)e 11 ( A 2 1 (t) 2 12 A 2 (t)sin() 11 A 1 (t)+ 2 12 A 2 2 (t)sin 2 () 2 11 ) dA 1 (t) = e 2 12 A 2 2 (t)sin 2 () 11 ∫ 1 0 A 1 (t)e 11 ( A 1 (t) 12 A 2 sin() 11 ) 2 dA 1 (t) (A.22) = e 2 12 A 2 2 (t)sin 2 () 11 ∫ 1 0 A 1 (t)e 11 (A 1 (t)) 2 dA 1 (t) (A.23) where≜ ( 12 A 2 (t)sin() 11 ) . Then, B(;A 2 (t); 11 ; 12 ) = e 2 12 A 2 2 (t)sin 2 () 11 ∫ 1 0 A 1 (t)e 1 2 (A 1 (t)) 2 1=2 11 dA 1 (t) (A.24) = e 2 12 r 2 2 sin 2 () 1 ∫ 1 (x+)e 1 2 x 2 1=2 1 dx (A.25) 122 wherex≜r. B(;A 2 (t); 11 ; 12 ) = e 2 12 A 2 2 (t)sin 2 () 11 [∫ 1 xe 1 2 x 2 1=2 11 dx+ ∫ 1 e 1 2 x 2 1=2 11 dx ] | {z } 1 2 11 e 11 2 + p 11 () (A.26) = e 2 12 A 2 2 (t)sin 2 () 11 [ 1 2 11 e 2 12 A 2 2 (t)sin 2 () 11 + 12 A 2 (t)sin() p 11 p 11 ( 12 A 2 (t)sin() 11 )] = 1 2 11 +e 2 12 A 2 2 (t)sin 2 () 11 12 A 2 (t)sin() p 11 p 11 ( 12 A 2 (t)sin() 11 ) where (:)isstandardnormalcumulativedistributionfunction. Hence,weget E [ e j∆ϕ(t) ] = 4 det(K z ) ∫ =2 =2 sin() ∫ 1 0 A 2 (t)e 22 A 2 2 (t) [ 1 2 11 +e 2 12 A 2 2 (t)sin 2 () 11 12 A 2 (t)sin() p 11 p 11 ( 12 A 2 (t)sin() 11 )] dA 2 (t)d = 4 12 p 11 det(K z ) ∫ 1 0 e 22 A 2 2 (t) ∫ =2 =2 e 2 12 A 2 2 (t)sin 2 () 11 2 12 A 2 2 (t)sin 2 () 11 ( 12 A 2 (t)sin() 11 ) d | {z } ≜C(A 2 (t); 11 ; 12 ) dA 2 (t) (A.27) 123 Bydefininga≜ 12 A 2 (t) 11 wesimplifyC(A 2 (t); 11 ; 12 ): C(A 2 (t); 11 ; 12 ) = 11 ∫ =2 =2 e 11a 2 sin 2 () (a 2 sin 2 ())(asin())d (A.28) = 11 a 2 [ ∫ 0 =2 e 11a 2 sin 2 () sin 2 ()(asin())d + ∫ =2 0 e 11a 2 sin 2 () sin 2 ()(asin())d ] = 11 a 2 ∫ 0 =2 e 11a 2 sin 2 () sin 2 ()(1(asin()))d + 11 a 2 ∫ =2 0 e 11a 2 sin 2 () sin 2 ()(asin())d = 11 a 2 ∫ 0 =2 e 11a 2 sin 2 () sin 2 ()d + 11 a 2 ∫ =2 0 e 11a 2 sin 2 () sin 2 () : 0 [(asin())+(asin())]d = 11 a 2 ∫ 0 =2 e 11a 2 sin 2 () sin 2 ()d (A.29) = 11 a 2 4 e 11 a 2 2 [ I 0 ( 11 a 2 2 ) +I 1 ( 11 a 2 2 )] (A.30) = 2 12 A 2 2 (t) 4 11 e 2 12 A 2 2 (t) 2 11 [ I 0 ( 2 12 A 2 2 (t) 2 11 ) +I 1 ( 2 12 A 2 2 (t) 2 11 )] (A.31) Hence, E [ e j∆ϕ(t) ] = p 12 det(K z ) 3=2 11 | {z } ≜D ∫ 1 0 A 2 2 (t)e 22 A 2 2 (t)+ 2 12 A 2 2 (t) 2 11 [ I 0 ( 2 12 A 2 2 (t) 2 11 ) +I 1 ( 2 12 A 2 2 (t) 2 11 )] dA 2 (t) = D ∫ 1 0 A 2 2 (t)e A 2 2 (t)d [ I 0 ( cA 2 2 (t) ) +I 1 ( cA 2 2 (t) )] dA 2 (t) (A.32) (A.33) whered≜ 22 2 12 2 11 andc≜ 2 12 2 11 . Then, E [ e j∆ϕ(t) ] = D [∫ 1 0 A 2 2 (t)e A 2 2 (t)d I 0 ( cA 2 2 (t) ) dA 2 (t)+ ∫ 1 0 A 2 (t) 2 e A 2 2 (t)d I 1 ( cA 2 2 (t) ) dA 2 (t) ] (A.34) 124 UsingMathematicatosolvethetwointegralsinEquA.34wefind: E [ e j∆ϕ(t) ] = D [ 2 F 1 ( 3 4 ; 5 4 ;1; c 2 d 2 ) p 4d 3=2 + 2 F 1 ( 5 4 ; 7 4 ;2; c 2 d 2 ) 3c p 16d 5=2 ] (A.35) = p 2 2 12 det(K z ) [ w 3=2 2 F 1 ( 3 4 ; 5 4 ;1;w 2 ) + 3 4 w 5=2 2 F 1 ( 5 4 ; 7 4 ;2;w 2 )] wherew = 2 12 2 11 22 2 12 and 2 F 1 representsahypergeometricfunction. Hence,weget E [ e j∆ϕ(t) ] = p 2 ( 1 11 22 2 12 )[ w 3=2 2 F 1 ( 3 4 ; 5 4 ;1;w 2 ) + 3 4 w 5=2 2 F 1 ( 5 4 ; 7 4 ;2;w 2 )] :(A.36) 125 AppendixB PartialPhaseSynchronization B.1 Conditional covariance matrix ofx conditioned on y FromBayesrule,wehave: p(xjy) = p(x;y) p(y) (B.1) = 1 M det(Kz) exp 8 < : [ x H y H ] K 1 z 2 4 x y 3 5 9 = ; 1 M2 det(Ky) exp { y H K 1 y y } : (B.2) Tosimplifythisratio,wedefineapartitionedformoftheinverseofthejointcovariance matrixsuchthat: 2 4 K x K xy K yx K y 3 5 | {z } Kz 2 4 A C H C B 3 5 | {z } K 1 z = 2 4 I 0 0 I 3 5 (B.3) whereIisMM identitymatrix. Thefollowingeaquationscanthenbeobtained: K x A+K xy C =I , K yx A =K y C: (B.4) 126 Now,theconditionaldistributioncanbewrittenas: p(xjy) = 1 M det(Kz) exp 8 < : [ x H y H ] 2 4 A C H C B 3 5 2 4 x y 3 5 9 = ; 1 M2 det(Ky) exp { y H K 1 y y } (B.5) = 1 m det(K z )det(K 1 x ) exp { [ x H Ax+2y H Cx+y H ( BK 1 y ) y ]} The exponential is a quadratic form in y and therefore conditional probability is also circular complex Gaussian. To determine the mean and covariance of the conditional distribution,weequate [ x H Ax+2y H Cx+y H ( BK 1 y ) y ] = (xm xjy ) H K 1 xjy (xm xjy ) (B.6) Then,wehave x H Ax = x H K 1 xjy x (B.7) y H Cx = x H K 1 xjy m xjy (B.8) y H ( BK 1 x ) y = y H K 1 xjy y (B.9) whicharevalidforallxandy implying K 1 xjy =A: (B.10) UsingtheequationsinB.4,weget: A 1 =K x K xy K 1 y K yx : (B.11) 127 Hence,theconditionalcovariancematrixcanbewrittenas K xjy =K x K xy K 1 y K yx : (B.12) Basically, in this setup the inverse of the conditional covariance is the upper left 2 2 sub-matrixoftheinversecovariancematrixwhichis 2 4 11 12 e j 12 21 e j 21 22 3 5 . 128 AppendixC ConnectivityinPresenceofLinear Mixing C.1 Linear Mixing Invariant Property of Lagged Coherenceinpresenceoftwosources Saytheobservedsignalsx m andx n inmodel4.3arerepresentedas: x m = a T m s (C.1) x n = a T n s (C.2) wherea m ≜ [ a mm a nm ] T ,a n ≜ [ a mn a nn ] T ands≜ [ s m s n ] T . Thecovari- ancestructureofthesourcevectorsisdefinedas: E [ ss H ] ≜R+jJ (C.3) where R ≜ R { E [ ss H ]} = 2 4 R mm R mn R nm R nn 3 5 (C.4) J ≜ I { E [ ss H ]} = 2 4 0 J mn J mn 0 3 5 =J mn 2 4 0 1 1 0 3 5 =J mn P: (C.5) 129 whereP≜ 2 4 0 1 1 0 3 5 . Using the definition of the LC given in equation 4.7, the LC betweens m ands n canbecomputedas: LC (s m ;s n ) = J mn √ R mm R nn R 2 mn : (C.6) AndtheLCbetweentheobservedsignalscanbecomputedas: LC (x m ;x n ) = IfE[x m x n ]g √ E[jx m j 2 ]E[jx n j 2 ]RfE[x m x n ]g 2 (C.7) = a T m Ja n √ a T m Ra m a T n Ra n (a T m Ra n ) 2 (C.8) = a T m Ja n √ a T m R[a m a T n a n a T m ]Ra n (C.9) = J mn a T m Pa n √ a T m R[a m a T n a n a T m ]R)a n (C.10) Notethat, [ a m a T n a n a T m ] = ( a T m Pa n ) P. Thenwehave: LC (x m ;x n ) = J mn a T m Pa n √ (a T m Pa n )(a T m RPRa n ) (C.11) = J mn a T m P pq a n √ (a T m Pa n ) 2 (R mm R nn R 2 mn ) (C.12) = J mn √ R mm R nn R 2 mn (C.13) whichissimplyequivalentto LC (s m ;s n ). 130 C.2 Proof of the relationship between Orthogonal CoherenceandtheStandardCoherence Thepartofthecomplexsignalthatcanbeinstantenouslyandlinearlybepredictedfrom x m ,i. e. x n∥m isdefinedas: x n∥m =R { E[x m x n ] E[x n x n ] } x m : (C.14) Thesignalx n orthogonalizedtox m ,i.e. x n?m ,canbederivedbysubtractingtheparallel signalcomponent: x n?m = x n x n∥m (C.15) = x n R { E[x m x n ] E[x n x n ] } x m : (C.16) Notethatthisisnotastandardprojection. Thestandardprojectionofx n ontox m where thetotalcross-coherenceiszerocanbedefinedas x ? n = E[x m x n ] E[x m x m ] x m (C.17) = R { E[x m x n ] E[x m x m ] } x m +jI { E[x m x n ] E[x m x m ] } x m (C.18) wherethecross-coherencebetweenx m andx ? n iszero: E [ x m (x n x ? n ) ] = E[x m x n ]E [ x m x ? n ] = 0: (C.19) 131 Bypluggingx ? n fromeq. C.18tothisequation,weget: E[x m x n ] =R { E[x m x n ] E[x m x m ] } E[x m x m ]+jI { E[x m x n ] E[x m x m ] } E[x m x m ] E[x m x n ]R { E[x m x n ] E[x m x m ] } E[x m x m ] | {z } E[xmx n?m ] =jI { E[x m x n ] E[x m x m ] } E[x m x m ] (C.20) implyingthatthecrossspectrumbetweenx m andx n?m is: E[x m x n?m ] =jIfE[x m x n ]g: (C.21) Inordertocomputecoherence,wealsoneedtocomputeauto-spectrumofx n?m : E[x n?m x n?m ] = E [ ( x n R { E[x m x n ] E[x n x n ] } x m ) 2 ] = E[x n x n ] RfE[x m x n ]g 2 E[x m x m ] (C.22) Hencethecoherencebetweenx m andx n?m canbesimplifiedto: C (x m ;x n?m ) = jIfE[x m x n ]g √ E[jx m j 2 ]E[jx n j 2 ]RfE[x m x n ]g 2 : (C.23) whichisidenticaltothelaggedcoherence. 132
Abstract (if available)
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Signal processing methods for interaction analysis of functional brain imaging data
PDF
Functional connectivity analysis and network identification in the human brain
PDF
Signal processing methods for brain connectivity
PDF
Applications of graph theory to brain connectivity analysis
PDF
Brain connectivity in epilepsy
PDF
Correction, coregistration and connectivity analysis of multi-contrast brain MRI
PDF
New theory and methods for accelerated MRI reconstruction
PDF
Shift-invariant autoregressive reconstruction for MRI
PDF
Multivariate statistical analysis of magnetoencephalography data
PDF
Functional real-time MRI of the upper airway
PDF
Improved brain dynamic contrast enhanced MRI using model-based reconstruction
PDF
Fast flexible dynamic three-dimensional magnetic resonance imaging
PDF
Optimization methods and algorithms for constrained magnetic resonance imaging
PDF
Diffusion MRI of the human brain: signal modeling and quantitative analysis
PDF
Causality and consistency in electrophysiological signals
PDF
High-dimensional magnetic resonance imaging of microstructure
PDF
Novel theoretical characterization and optimization of experimental efficiency for diffusion MRI
PDF
Geometric methods of image registration and analysis
PDF
Improving the sensitivity and spatial coverage of cardiac arterial spin labeling for assessment of coronary artery disease
PDF
Theory and simulation of diffusion magnetic resonance imaging on brain's white matter
Asset Metadata
Creator
Aydore, Sergul
(author)
Core Title
Measuring functional connectivity of the brain
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Electrical Engineering
Publication Date
03/04/2015
Defense Date
08/14/2014
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
brain,EEG,functional connectivity,MEG,OAI-PMH Harvest
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Leahy, Richard M. (
committee chair
), Haldar, Justin P. (
committee member
), Nayak, Krishna S. (
committee member
)
Creator Email
sergul@usc.edu,sergulaydore@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c3-468390
Unique identifier
UC11286552
Identifier
etd-AydoreSerg-2881.pdf (filename),usctheses-c3-468390 (legacy record id)
Legacy Identifier
etd-AydoreSerg-2881.pdf
Dmrecord
468390
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Aydore, Sergul
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
brain
EEG
functional connectivity
MEG