Augmenting interictal mapping with neurovascular coupling biomarkers by structured factorization of epileptic EEG and fMRI data

Delft University of Technology

Augmenting interictal mapping with neurovascular coupling biomarkers by structured

factorization of epileptic EEG and fMRI data

Van Eyndhoven, Simon ; Dupont, Patrick; Tousseyn, Simon ; Vervliet, Nico; Van Paesschen , Wim; Van

Huffel, Sabine; Hunyadi, Borbala

DOI

10.1016/j.neuroimage.2020.117652

Publication date

2021

Document Version

Final published version

Published in

NeuroImage

Citation (APA)

Van Eyndhoven, S., Dupont, P., Tousseyn, S., Vervliet, N., Van Paesschen , W., Van Huffel, S., & Hunyadi,

B. (2021). Augmenting interictal mapping with neurovascular coupling biomarkers by structured factorization

of epileptic EEG and fMRI data. NeuroImage, 228, 1-20. [117652].

https://doi.org/10.1016/j.neuroimage.2020.117652

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

NeuroImage228(2021)117652

ContentslistsavailableatScienceDirect

NeuroImage

journalhomepage:www.elsevier.com/locate/neuroimage

Augmenting

interictal

mapping

with

neurovascular

coupling

biomarkers

by

structured

factorization

of

epileptic

EEG

and

fMRI

data

Simon

Van

Eyndhoven

_,

_Patrick

_Dupont

_,

_Simon

_Tousseyn

_,

_Nico

_Vervliet

_,

_Wim

_Van

Paesschen

_,

_Sabine

_Van

_Huffel

_,

_Borbála

_Hunyadi

a Department of Electrical Engineering (ESAT), STADIUS Center for Dynamical Systems, Signal Processing and Data Analytics, KU Leuven, Belgium b Laboratory for Cognitive Neurology, Department of Neurosciences, KU Leuven, Leuven, Belgium

c Leuven Brain Institute, Leuven, Belgium

d Academic Center for Epileptology, Kempenhaeghe and Maastricht UMC+, Heeze, the Netherlands e Laboratory for Epilepsy Research, KU Leuven, Leuven, Belgium

f Department of Neurology, University Hospitals Leuven, Leuven, Belgium

g Circuits and Systems Group (CAS), Department of Microelectronics, Delft University of Technology, Delft, the Netherlands

a

r

t

i

c

l

e

i

n

f

o

Keywords:

EEG-fMRI

Blind source separation Tensor factorization Interictal epileptic discharge Neurovascular coupling Hemodynamic response function

a

b

s

t

r

a

c

t

EEG-correlatedfMRIanalysisiswidelyusedtodetectregionalBOLDfluctuationsthataresynchronizedto inter-ictalepilepticdischarges,whichcanprovideevidenceforlocalizingtheictalonsetzone.However,thetypical, asymmetricalandmass-univariateapproachcannotcapturetheinherent,higherorderstructureintheEEGdata, normultivariaterelationsinthefMRIdata,anditisnontrivialtoaccuratelyhandlevaryingneurovascular cou-plingoverpatientsandbrainregions.Weaimtoovercomethesedrawbacksinadata-drivenmannerbymeans ofanovelstructuredmatrix-tensorfactorization:thesingle-subjectEEGdata(representedasathird-order spec-trogramtensor)andfMRIdata(representedasaspatiotemporalBOLDsignalmatrix)arejointlydecomposed intoasuperpositionofseveralsources,characterizedbyspace-time-frequencyprofiles.Inthesharedtemporal mode,Toeplitz-structuredfactorsaccountforaspatiallyspecific,neurovascular‘bridge’betweentheEEGand fMRItemporalfluctuations,capturingthehemodynamicresponse’svariabilityoverbrainregions.Byanalyzing interictaldatafromtwelvepatients,weshowthattheextractedsourcesignaturesprovideasensitivelocalization oftheictalonsetzone(10/12).Moreover,complementarypartsoftheIOZcanbeuncoveredbyinspectingthose regionswiththemostdeviantneurovascularcoupling,asquantifiedbytwoentropy-likemetricsofthe hemo-dynamicresponsefunctionwaveforms(9/12).Hence,thismultivariate,multimodalfactorizationprovidestwo usefulsetsofEEG-fMRIbiomarkers,whichcanassistthepresurgicalevaluationofepilepsy.Wemakeallcode requiredtoperformthecomputationsavailableathttps://github.com/svaneynd/structured-cmtf.

1. Introduction

Refractoryepilepsy is aneurologicaldisordersufferedby30%of approximately50 millionepilepsypatientsworldwide(WorldHealth Organization,2019),inwhichseizurescannotadequatelybecontrolled byanti-epilepticmedication.Inthepreparationoftreatmentvia resec-tivesurgery,interictalepilepticdischarges(IEDs)canbelocalizedinthe brainwithsimultaneousEEG–fMRI,whichprovidesagoodsurrogatefor mappingtheseizureonsetzone(Anetal.,2013;Grouilleretal.,2011; vanHoudtetal.,2013;Khooetal.,2017;Lemieuxetal.,2001; Thorn-tonetal.,2010;Vulliemozetal.,2009;Zijlmansetal.,2007).Thisis oftenconductedviaEEG-correlatedfMRIanalysis,whereinareference temporalrepresentationoftheIEDsisusedtointerrogateallbrain re-gions’bloodoxygenleveldependent(BOLD)signalsforsignificant

cor-∗_{Correspondingauthor.}

E-mailaddresses:simon.vaneyndhoven@gmail.com,simon.vaneyndhoven@gmail.com(S.VanEyndhoven).

relations;voxelsforwhichastatisticalthresholdisexceededcanthen beconsideredpartoftheepilepticbrainnetwork,alongwhichepileptic seizuresaregeneratedandpropagated(Gotman,2008;Lemieuxetal., 2001;Salek-Haddadietal.,2003;Thorntonetal.,2010;Zijlmansetal., 2007).

Since its inception, theworkhorsefor conducting EEG-correlated fMRI analysis hasbeen the general linear model(GLM) framework (Fristonetal.,1994;PolineandBrett,2012;Salek-Haddadietal.,2006). Overthepastyears,ithasbecomeclearthatusingtheGLMcomeswith severalhurdles,relatedtothemanymodelingassumptions,thatmay reduceitssensitivityorspecificity(increasingTypeIerrors)when vi-olated(Lindquistet al.,2009; Monti,2011; PolineandBrett, 2012). Remediesforseveraloftheseissuesarenotyetwidelyapplied,orare notyetavailable.

https://doi.org/10.1016/j.neuroimage.2020.117652

Received15April2020;Receivedinrevisedform28November2020;Accepted4December2020 Availableonline24December2020

Firstofall,theadoptionofarelevantrepresentationofIED occur-rencestoconstructaregressorforthedesignmatrixhasprovenvitalto thesensitivity.ThisaspecthasbeeninvestigatedinRosaetal.(2010),

Murtaetal.(2015),Abreuetal.(2018),VanEyndhovenetal.(2019a). Inpreviouswork(VanEyndhovenetal.,2019a),weaddressedthis is-suebypre-enhancingtheEEGsignalsusingaspatiotemporalfilterthatis tunedtomaximizethesignal-to-noiseratio(SNR)ofIEDswithrespect tothebackgroundEEG.Wehaveshownthattakingthetime-varying powerofthefilteredEEGleadstoarobustregressor,whichismore per-formantthanmanyothertypesofregressors,includingthosebasedon stickfunctions(Lemieuxetal.,2001;Salek-Haddadietal.,2006),ICA (Abreu etal.,2016; Formaggioet al.,2011) orEEGsynchronization (Abreuetal.,2018).

Modelmismatchmayoccurduetotheunknownneurovascular cou-plingfrom electrophysiologicalphenomenameasured onthe EEGto hemodynamicvariationscapturedbytheBOLDsignals.Inmanypapers onEEG-correlatedfMRI, acanonicalhemodynamicresponsefunction (HRF)basedontwogammadensityfunctionsisusedtotranslate IED-relatedtemporaldynamicstoBOLDfluctuations(Fristonetal.,1998). However,thereisinsurmountableevidencethattheHRFisnotfixed, butvariessubstantiallyoversubjects(Aguirreetal.,1998),overbrain regions(Handwerkeretal.,2004),withage(Jacobsetal.,2008),or evenwithstresslevel(Elbauetal.,2018).Forthediseasedbrain,this issuemaybe evengreater:i.e.,additionalvariation,e.g.in brain ar-easinvolvedintheepilepticnetwork,hasbeenobservedcomparedto healthycontrols(Bénaretal.,2002;Grouilleretal.,2010;vanHoudt etal.,2013;Jacobsetal.,2009;Lemieuxetal.,2008).Plentyofprevious researchhasshownthatfailingtoaccountforthisvariabilitymaylead tosubstantialbiasandincreasedvarianceoftheestimatedactivation, whichinturninflatesTypeIand/orTypeIIerrorrates(Calhounetal., 2004;Lindquistetal.,2009;LindquistandWager,2007;Monti,2011). Severalmethodshavebeendevisedtodealwiththisvariability.A widelyusedapproachistomodeltheHRFasalinearcombinationof severalbasisfunctions.Somepopularchoicesforthesebases,whichare alsosupportedbyopensourcetoolboxeslikeSPMarethe‘informed ba-sisset’(Fristonetal.,1998),consistingoftheHRFplusitsderivative w.r.t.timeanditsderivativew.r.t.thedispersionparameter(leading toaTaylor-likeextensionwhichcancaptureslightchangesinpeak on-setandwidth),andthefiniteimpulsereponse(FIR)basisset,inwhich everybasisfunctionfitsexactlyonesampleoftheHRFineveryvoxel (Aguirreetal.,1998;Glover,1999).Otherresearchershaveaimedto findabasissetbycomputingalow-dimensionalsubspaceofalargeset of‘reasonable’HRFs(Woolrichetal.,2004)orbyfittingnonlinear func-tionstogivenfMRIdata(LindquistandWager,2007;VanEyndhoven etal.,2017).Alternatively,multiplecopiesofastandardHRF,which differonlyintheirpeaklatencies,canbeused(Bagshawetal.,2004). Finally,approachesexistthataimtobeimmunetodifferencesin neu-rovascularcoupling,suchasthosebasedonmutualinformation(MI), whichdoesnotrelyonanypredefinedmodelorevenlinearityofthe HRF(Caballero-Gaudesetal.,2013;OstwaldandBagshaw,2011). Per-hapssurprisingly,theauthorsofCaballero-Gaudesetal.(2013)found thattheresultsbasedonMIwereoftenverysimilartothosebasedon theinformedbasisset,leadingtotheconclusionthattheassumptionof alineartime-invariantsystem,asdescribedbytheconvolutionwithan appropriateHRF,issufficientlyaccurate.Instead,itmaybeusefulto notmakeabstractionofthevariableneurovascularcoupling,butrather consideritasanadditionalbiomarkertolocalizeepileptogeniczones (vanHoudtetal.,2013).Indeed,inseveralstudies,HRFsthatdeviate fromthecanonicalmodelwerefoundinregionsoftheepilepticnetwork (Bénaretal.,2002;Hawcoetal.,2007;vanHoudtetal.,2013;Jacobs etal.,2009;Lemieuxetal.,2008;Moeller etal.,2008;Pittauetal., 2011).Severalhypotheseshavebeenpostulatedtoexplainthis varabil-ity,includingalteredautoregulationduetohighermetabolicdemand following(inter)ictalevents(Schwartz,2007),vascularreorganization neartheepileptogenicregion(Rigauetal.,2007),ortheexistenceof pre-spikechangesinneuro-electricalactivitywhicharenotvisibleon

EEGandwhichculminateintheIED(Jacobsetal.,2009).Itisthusan opportunitytomapnotonlyregionswithstatisticallysignificantBOLD changesinresponsetoIEDs,butalsothespatialmodulationoftheHRF waveformitself, inorder todiscoverregions whereanaffectedHRF shapemayprovideadditionalevidencetowardstheepilepticonset.

Thepreviousconsiderationsindicatethatitisdifficulttomeetall assumptionsinthegenerallinearmodel,whichmaycompromise infer-ence power(Handwerkeretal., 2004; Lindquistetal., 2009; Monti, 2011). Data-driven alternatives may relieve this burden, since they adapttothecomplexity ofthedatamoreeasilycomparedto model-based approaches,andare especiallysuitedfor exploratoryanalyses (Mantinietal.,2007;Marečeketal.,2016).BlindSourceSeparation (BSS)techniquesconsiderEEGand/orfMRIdatatobeasuperposition of several‘sources’ofphysiological activityandnonphysiological in-fluences. Basedontheobserveddataalone,BSStechniquesareused to estimate boththe sourcesandthe mixingsystem, bymeansof a factorization ofthedatainto two(ormore)factormatrices,holding sourcesormixingprofilesalongthecolumns. Theynaturallyallowa symmetricaltreatmentofEEGandfMRIdata,enablingtruefusionof bothmodalities(Calhounetal.,2009;Lahatetal.,2015;Valdes-Sosa etal.,2009),whichisincontrasttoEEG-correlatedfMRI,where EEG-derivedIEDsinformthefMRIanalysis.Whiletheinformation-theoretic approach in Caballero-Gaudes et al.(2013) alsoshares this symme-tryfeature,itpurposelyavoidstheestimationofHRFs, whichisour goalhere.Furthermore,BSStechniquesnaturallyaccommodate higher-order representationsof thedatain theformof tensorsormultiway arrays,whichcancapturetherichstructureinthedata.Indeed, mea-surementsofbrainactivityinherentlyvaryalongseveralmodes (sub-jects,EEGchannels,frequency,time,...),whichcannotberepresented usingmatrix-basedtechniqueslikeICAwithoutlossofstructureor in-formation (Acar etal., 2007; Lahat etal., 2015; Sidiropoulos etal., 2017).Tensor-basedBSStechniqueshavebeenusedtomineunimodal EEGdatabydecomposingthird-order spectrograms(channels×time points×waveletscales)intoseveral‘atoms’(alsocoined‘components’ or ‘sources’),each withadistinct spatial,temporalandspectral pro-file/signature(Marečeketal.,2016; Miwakeichietal., 2004;Mørup etal.,2006),withsuccessfulapplicationinseizureEEGanalysis(Acar etal.,2007;DeVosetal.,2007).WhileatensorextensionofICAfor groupfMRIdata(intheformof subjects×timepoints ×voxels) ex-ists (BeckmannandSmith,2005),matrixrepresentations offMRI re-main dominantforsingle-subject analyses.Moreover, sucha tensor-based extension implicitlyassumes that sourceshave thesame time course oversubjects,which isnotanadequatemodelforIED occur-rences,norresting-statefluctuations.CoupledBSStechniquescan esti-matecomponentswhicharesharedbetweenbothmodalities,providing acharacterizationinbothdomains(Hunyadietal.,2017).For exam-ple, inAcar etal.(2017), Acaretal.(2019), Hunyadietal.(2016),

Chatzichristosetal.(2018),multi-subjectEEGandfMRIdatahavebeen analyzedusingcoupledmatrix-tensorfactorization(CMTF),whereinthe ‘subjects’factorissharedbetweentheEEGtrilineartensor decompo-sition andthefMRImatrixdecomposition.InHunyadi etal.(2016), theresultingfactorsignaturesrevealedonsetandpropagationzonesof aninterictalepilepticnetworkthatwascommonoverpatients,aswell asthemodulationofthedefault-modenetwork(DMN)activity.Also single-subjectdatacanbedecomposedintodistinctcomponents,using a sharedtemporalfactorforEEG andfMRI.This requirestheuse of amodeloftheneurovascularcoupling,toensuretemporalalignment ofEEGandBOLDdynamics.InMartínez-Montesetal.(2004),afixed canonicalHRFwasused,followedbymultiwaypartialleastsquaresto extractcomponentswithspatial,temporal,andspectralsignatures.In previous work,weproposed anextensiontothistechnique, wherea subject-specificHRFisco-estimatedfromtheavailabledata,alongwith thecomponents(VanEyndhovenetal.,2017).

Inthispaper,weextendthislattertechniqueinordertoaccountnot onlyforsubject-wisevariationoftheHRF,butalsocapturevariations overbrainregions.ThisresultsinahighlystructuredCMTF(sCMTF)of

S. Van Eyndhoven, P. Dupont, S. Tousseyn et al. NeuroImage 228 (2021) 117652 Table1

Clinicalpatientdata.

patient gender ictal onset zone etiology surgery

ILAE outcome

follow-up

time (y) # IEDs

# TRs (# sessions)

# EEG

channels IED loc. p01 F L temporal HS temporal lobe

resection 3 5 15 540 (1) 29 F7–T1 p02 F L parietal FCD partial lesionectomy 4 5 663 1620 (3) 29 Pz p03 F R parieto-occipito- temporal Sturge-Weber 105 1080 (4) 21 F8 p04 M R temporal unknown 825 1620 (3) 21 F8–T4

p05 F L anterior temporal HS temporal lobe

resection 1 8 117 1080 (3) 29 F7–T1

p06 F R frontal FCD partial

lesionectomy

5 2 640 1080 (3) 29 Cz–C4

p07 F L anterior temporal DNET temporal lobe resection

1 4 126 1080 (4) 29 F7–T1

p08 M L temporo-parietal unknown overlap eloquent cx 11 1080 (4) 21 T5 p09 F L temporo-occipital FCD overlap eloquent cx 1815 1620 (3) 29 T3–T5 p10 F R temporal HS refused 226 540 (1) 29 F8–T2

p11 M L anterior temporal HS temporal lobe resection

1 6 6 1080 (2) 29 F7–T1

p12 F R temporal CNS infection refused 966 1350 (5) 27 T4

Abbreviations:F=female,M=male,L=left,R=right,CNS=centralnervoussystem,DNET=dysembryoplasticneuroepihelialtumor,FCD=focalcortical dysplasia,HS=hippocampalsclerosis,cx=cortex,IED=interictalepilepticdischarge,TR=repetitiontime,IEDloc.=localizationoftheIEDonEEG.

theinterictalmultimodaldata,inwhichHRFbasisfunctionsand spa-tialweightingcoefficientsareestimatedalongwithspatial,spectraland temporalsignaturesofcomponents.BypreprocessingtheEEGusingthe data-drivenfiltersfromVanEyndhovenetal.(2019a),weaimto max-imizethesensitivityinmappingtheinterictaldischarges.Weanalyze whethertheestimatedspatialmodulationoftheHRFwaveformisa vi-ablebiomarkerwhenlocalizingtheictalonsetzone,besidestheBOLD spatialsignaturesthemselves.

2. Methodsandmaterials 2.1. Patientgroup

We use data of twelve patients, whom we previously studiedin

Tousseynetal.(2014a),Tousseynetal.(2014b),Tousseynetal.(2015),

Hunyadietal.(2015),VanEyndhovenetal.(2019a). Thesepatients wereselectedbasedonthefollowingcriteria:(1)theywereadultswhich underwent presurgicalevaluation for refractory focalepilepsy using EEG–fMRI,andforwhichtherewasconcordanceofalltheavailable clin-icalevidenceregardingtheepilepticfocus;(2)subtractionictal single-photonemission tomography(SPECT)coregistered toMRI (SISCOM) imageswereavailableforallpatients,aswellaspost-surgeryMRIscans whenpatientswereseizure-free(internationalleagueagainstepilepsy (ILAE)outcome classification1-3 (1,completelyseizure-free;2, only auras;3,onetothreeseizuredaysperyear±auras;4,fourseizuredays peryearto50%reductionofbaselineseizuredays±auras;5,<50% reductionofbaselineseizuredaysto100%increaseofbaselineseizure days±auras;6,morethan100%increaseofbaselineseizuredays± auras));(3)IEDswererecordedduringtheEEG–fMRIrecordingsession. Thisstudywascarriedoutinaccordancewiththerecommendations oftheInternationalConferenceonHarmonizationguidelinesonGood ClinicalPracticewithwritten informedconsentfromallsubjects. All subjectsgavewritteninformedconsentinaccordancewiththe Decla-rationofHelsinki,fortheirdatatobeused inthisstudy, butnotto bemadepubliclyavailable.Theprotocolwasapprovedbythe Medi-calEthicsCommitteeoftheUniversityHospitalsKULeuven.Forthe patients’completeclinicaldata,werefertoTable1.

2.2. Dataacquisitionandpreprocessing

Functional MRI data were acquired on one of two 3TMR scan-ners(AchievaTXwitha32-channelheadcoilandInteraAchievawith aneight-channelheadcoil,PhilipsMedicalSystems,Best,The Nether-lands)withanechotime(TE)of33ms,arepetitiontime(TR)of ei-ther 2.2or2.5s, andavoxelsizeof 2.6×3×2.6mm3_{.EEG data}

wererecordedaccordingtotheinternational10–20systemusing MR-compatiblecaps,sampledat5kHz,withCzreference.TheEEGsignals wereband-passfilteredofflinebetween1-50Hz,gradientartifactswere removedusingtheBergenplug-in(BergenfMRIGroup,Bergen,Norway) forEEGLAB(Moosmannetal.,2009)andpulseartifactsweresubtracted withtheBrainVisionAnalyzersoftware(BrainProducts,Munich, Ger-many)(Allenetal.,1998).Thesignalofeverychannelwasdividedby itsstandarddeviation.Twoneurologistssubsequentlyinspectedand an-notatedtheEEGsignalsforIEDs.

ThefMRIimageswererealigned,slice-timecorrectedand normal-izedtoMNIspace,resampledtoavoxelsizeof2×2×2mm3_,and

smoothedusingaGaussiankernelof6mmfullwidthathalfmaximum (FWHM).These processingstepswerecarriedoutusingSPM8 (Func-tionalImaging Laboratory, Wellcome Centerfor Human Neuroimag-ing,UniversityCollegeLondon,UK)(Fristonetal.,1994).Wereferthe readertoTousseynetal.(2014a)foradetaileddescriptionofthese pre-processingsteps.

WeregressoutcovariatesofnointerestfromthefMRIdata.These include:thesixmotion-correction parameters(plustheirsquares and derivatives);boxcarregressorsatmomentsofsubstantialscan-to-scan headmovement(largerthan1mmbasedonthetranslation parame-ters);thefirstfiveprincipalcomponentsoftheBOLDtimeserieswithin thecerebrospinalfluidandwhitematterregions(Behzadietal.,2007). Subsequently,theBOLDtimeseriesarefilteredbetween0.008–0.20Hz usingtheCONNtoolbox(Whitfield-GabrieliandNieto-Castanon,2012). Forananalysisoftheeffectoftheorderingofthesepreprocessingsteps, werefertothesupplementarymaterial.

The dimensionalityof thefMRI data is reduced by meansof an anatomicalparcellationofthebrain.Theinitial79× 95× 68imagesare segmentedintoregions-of-interest(ROIs)accordingtotheBrainnetome

atlas,whichconsistsof246parcelsinthegreymatter(Fanetal.,2016). ForeveryROI,oneBOLDtimeseriesisconstructedastheaverageof thetimeseriesofallvoxelswithintheROI.Ifmultipleacquisitionruns (withinthesamerecordingsession)hadbeendone,theEEGandfMRI dataof the different runs aretemporally concatenated. Further cus-tomizedpreprocessingstepsaretreatedinSections2.3and2.4.

2.3. Multi-channelWienerfilteringforspatio-temporalEEGenhancement

Inpreviouswork(VanEyndhovenetal.,2019a),wehaveshownthat pre-enhancingtheEEGsignalsusingadata-driven,spatiotemporalfilter thatistunedtomaximizethesignal-to-noiseratio(SNR)ofIEDswith re-specttothebackgroundEEGandartifacts,leadstoaBOLDpredictorthat ismoreperformantthanmanyotherpredictors,includingthosebased onsimplestickfunctions(Lemieuxetal.,2001;Salek-Haddadietal., 2006),ICA(Abreuetal.,2016;Formaggioetal.,2011) orEEG syn-chronization(Abreuetal.,2018).Thispre-enhancementstrategybased onmulti-channelWienerfilters(MWF)haserror-correctingcapabilities andproducesanIEDrepresentationthatimprovesthelocalization sen-sitivityofEEG-correlatedfMRI(VanEyndhovenetal.,2019a).

Inbrief,theMWFisestimatedbyfirstperformingtime-delay embed-dingofthemulti-channelEEGsignals𝐱[𝑡]∈ℝ𝐼𝑚_,leadingtoanextended multi-channel,multi-lagsignal̃𝐱[𝑡]∈ℝ2𝐼𝑚𝜏+𝐼𝑚as

̃𝐱= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝐱[𝑡−𝜏] ⋮ 𝐱[𝑡] ⋮ 𝐱[𝑡+𝜏] ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (1)

andsubsequentlycomputingthefiltercoefficientsas

̂ 𝐖=𝐑−1

𝑥𝑥(𝐑𝑥𝑥−𝐑𝑛𝑛), (2)

where 𝐑𝑥𝑥=E{̃𝐱̃𝐱T|𝐻=1} is the covariance matrix of the

EEG observed during annotated IED segments (𝐻=1), and

𝐑𝑛𝑛=E{̃𝐱̃𝐱T|𝐻=0} is the covariance matrix of the EEG outside

ofIEDsegments(𝐻=0).Forthefullderivation,wereferthereaderto (Somersetal.,2018;VanEyndhovenetal.,2019a).TheEEGsignals arethenfilteredas𝐖̂T_{̃𝐱. Duetotheextensionwithlaggedcopiesof}

thesignals,channel-specificfiniteimpulse responsefiltersarefound. Hence,𝐖̂T_{̃𝐱isasetofspatiotemporallyfilteredoutputsignals,inwhich}

IED-like waveformsarepreserved whileother waveforms,whichare notspecifictoepilepsy,aresupressed1_.

WetraintheMWFforeachpatientindividually, afterembedding theEEGsignalsusing𝜏 =4positiveandnegativelags2_,andcompute thefinalfilterusingthegeneralizedeigenvaluedecomposition,which ensuresthepositivedefinitenesspropertyofthesubtractedcovariance matrixin(2)(Somersetal.,2018).

2.4. Higher-orderdatarepresentation

Topreservetheintrinsicmultiwaynatureofthedata,werepresent thepreprocessedEEG andfMRI asatensorandmatrixrespectively, whicharesubsequentlyfactorizedjointly.Thisapproachdiffers from themass-univariatetreatmentinthetraditionalGLM,whereeachvoxel istreatedindividually,andonly‘flattened’EEGtimecoursescanbe en-teredasregressors.Sinceepilepsyismanifestedwithconsiderable vari-abilitybetweenpatients,wehandlethemultimodaldataofeachpatient separately.

1 _{Toretrievefilteredversionsoftheoriginalsetofchannelsonly(andnotof}

theirtime-delayembeddedcopies,whichwouldberedundant),onlythemiddle

𝐼𝑚columnsof𝐖̂ areused(cfr.(1)).

2 _{WeobservedinVanEyndhovenetal.(2019a)thatformostpatients,the}

outputSNRsaturatesaroundthisvalue,correspondingtoanintervalof-16ms to+16ms.

2.4.1. Spatio-temporal-spectraltensorrepresentationofEEG

Weadopt atensorizationstrategybased ontime-frequency trans-formationoftheEEGdatatothird-orderspectrograms(timepoints× frequencies×channels).Afterthepre-enhancementstepdescribedin

Section 2.3, we create aspectrogram using theThomson multitaper method,appliedonnonoverlappingEEGsegmentswithalengthequal toonerepetitiontime(TR)ofthefMRIacquisition.ThesquaredFourier magnitudesareaveragedinto1Hzbins,from1Hzto40Hz.Hence,for everyEEGchannel,weobtainaspectrogramwhichissynchronizedto thefMRItimeseries.Thetimepoints×frequencies×channels spectro-gram∈ℝ𝐼𝑠×𝐼_𝑔×𝐼_{𝑚isfurthernormalizedasdescribedinA.1,to}

equal-izetheinfluenceofeachchannelandeachfrequency,andtofocuson relativesignalincreasesordecreases(Marečeketal.,2017;2016)

2.4.2. Spatio-temporalmatrixrepresentationoffMRI

TheaverageBOLDtimeseriesarestackedinatimepoints×ROIs matrix𝐘∈ℝ𝐼𝑠×𝐼𝑣_,where𝐼𝑣=246ROIs.WenormalizeeachROI’stime

seriesbysubtractingitsmeananddividingbyitsstandarddeviation.

2.4.3. Neurovascularcouplinginthetemporalmode

EEGandfMRIdataareacquiredsimultaneouslypersubject,andare thusnaturallycoregisteredalongthe‘time’mode.Thisiscapturedina temporalfactormatrixthatiscommonbetweentheEEGfactorization andthefMRIfactorization.However,theelectrophysiologicalchanges thatarepickedupbyEEGvaryonamuchmorerapidtimescalethan thesluggishBOLDfluctuationsthat(indirectly)correspondtothesame neuralprocess.Theneurovascularcouplingthatdescribestherelation betweenthesetwocomplementarysignalscanbedescribedbya convo-lutionwithanHRF3_.

Inpreviouswork,wedevelopedaCMTFmodelinwhichtheHRF itselfisparametricallyestimatedfromthedata(VanEyndhovenetal., 2017),andamatrixmultiplicationwithToeplitzstructureimplements theHRFconvolution,asproposedinValdes-Sosaetal.(2009).Inthe samepaper,wehintedtowardsanextensionbasedonmultiple basis functionstoaccountforthevariabilityoftheHRFoverbrainregions. Inthefollowing,weassumethatthetimecourseofeachEEGsource isconvolvedwithanaprioriunknown,ROI-specificHRF,whichisa superpositionof𝐾 parametrizedbasisfunctions,whichleadstoa mod-elledcontributionof thissourcetotheROI’sBOLDsignal. Hence,in everyROI𝑖𝑣,themodeled(unscaled)BOLDtimecourse𝐳𝑖(𝑣𝑟)ofthe𝑟-th

neuralsourceis 𝐳(𝑟) 𝑖𝑣 =𝐇𝑖𝑣𝐬𝑟 (𝑟=1…𝑅) (3) = 𝐾 ∑ 𝑘=1 𝑏𝑘,𝑖𝑣𝐇𝑘𝐬𝑟 (4) = 𝐾 ∑ 𝑘=1 𝑏𝑘,𝑖𝑣 ( 𝐡𝑘)𝐬𝑟 (5) = 𝐾 ∑ 𝑘=1 𝑏𝑘,𝑖𝑣 ( (𝜽𝑘₎)_𝐬 𝑟. (6)

Here,𝐬𝑟isafactorvectorholdingthetimecourseofthe𝑟-thEEGsource;

 isanoperatorthattransformsasetofparameters 𝜽(𝑘)_{intoafullHRF,}

representedasavector𝐡𝑘; isanoperatorthattransforms𝐡𝑘intoa

Toeplitzmatrix𝐇𝑘bypopulatingthemainandlowerdiagonalswith

theHRFsamples(seealsoA.2);𝑏_{𝑘,𝑖𝑣}istheweightforthe𝑘-thHRFbasis functioninthe𝑖𝑣-thROI;𝐇𝑖𝑣 istheToeplitzmatrixholdingthetotal

HRFinthe𝑖𝑣-thROI.ThisoperationisclarifiedinFig.1b.

3_{Inthispaper,weusetheterm‘neurovascularcoupling’todescribethe}

cou-plingcharacteristicbetweenEEGandfMRItemporaldynamics,andmakethe silentassumptionthatthischaracteristicisaproxy/surrogatefor ‘neurovascu-larcoupling’asitisunderstoodinneuroscience:themodelthatdescribesBOLD changesinresponsetoelectricalneural‘events’,whichtaketheformoflocal fieldpotentialsatthesynapses.

S. Van Eyndhoven, P. Dupont, S. Tousseyn et al. NeuroImage 228 (2021) 117652

Fig.1. Structuredcoupledmatrix-tensorfactorization(sCMTF)ofEEGandfMRIdatacanrevealneuralsourcesthatareencodedinbothmodalities,aswellascapture thevaryingneurovascularcouplingbetweentheelectrophysiologicalandBOLDchanges.(a)TheEEGsignalsvaryovertimepoints×frequencies×electrodes.The resultingthird-orderspectrogramtensorisfactorizedaccordingto(8)into𝑅rank-1components,whicheachconsistofatemporalsignature𝐬𝑟,aspectralsignature

𝐠𝑟andaspatialsignature𝐦𝑟.(b)ThefMRIdataconsistoftheaverageBOLDsignalindifferentbrainparcelsorregionsofinterest(ROIs),representedinatime

points×ROImatrix𝐘,andarefactorizedaccordingto(11).NeurovascularcouplingismodeledasaconvolutionoftheEEGtemporaldynamicswithaROI-specific hemodynamicresponsefunction(HRF),asin(11)–(13).Inthisexample,eachlocalHRFisrepresentedasalinearcombination(encodedbycoefficients𝐛𝑘)of

𝐾=3optimizedbasisfunctions,eachpopulatingaToeplitzmatrix𝐇𝑘whichimplementsaconvolutionthroughmatrixmultiplicationwiththetemporalsignatures

𝐬𝑟.Afterwards,eachsmoothedcomponent𝑟isspatiallyweightedbyasignature𝐯𝑟.Thisisaccomplishedbytheelementwiseproduct𝐛𝑘∗𝐯𝑟oftheHRFbasis

function-specificcoefficients𝐛𝑘andthecomponent-specificamplitudes𝐯𝑟.

Thistimecourse𝐳(𝑟)

𝑖𝑣 isconceptuallyequivalenttoaregressorinthe

GLM’sdesignmatrix.WetreattheHRFparametersets𝜽(𝑘)_{,𝑘=1…𝐾}

as unknown variables, which need tobe fitted tothe dataat hand (LindquistandWager,2007).Byparametrizingeachbasisfunction,we embedprotectionagainstnonsensicalHRFshapes,andagainst overfit-ting,sincethenumberofparameterstobeestimatedisgreatlyreduced comparedtotheFIRbasisinGlover(1999),Aguirreetal.(1998).We employa double-gammaHRF,i.e.,each HRFbasis function𝑘 is de-scribedby fiveparameters asℎ𝑘(𝑡)=𝑓(𝑡;𝜽)=Γ(𝜃1)−1⋅ 𝜃2𝜃1𝑡𝜃1

−1_𝑒−𝜃2𝑡₋

𝜃5Γ(𝜃3)−1⋅ 𝜃𝜃43𝑡𝜃3

−1_𝑒−𝜃4𝑡,_{whereweomitthesuperscript}(𝑘)_fromthe

pa-rameters𝜽 tonotoverloadthenotation.

2.5. Coupledmatrix-tensorfactorizationofEEGandfMRI

Aftertensorization,wejointlydecomposetheEEGtensorandthe fMRImatrix𝐘intoasetofdistinctsources.

Thethird-orderEEG spectrogramis approximatedbya sumof 𝑅

rank-1 termsaccordingtothetrilinearcanonicalpolyadic decompo-sition(CPD)(also referredtoasParallelFactorAnalysis(PARAFAC)) as in Miwakeichi et al. (2004), Mareček et al. (2016), Martínez-Montesetal.(2004),VanEyndhovenetal.(2017).Eachrank-1term

𝐬𝑟◦ 𝐠𝑟◦ 𝐦𝑟describesasource(also called‘component’)intermsofan

outerproduct(◦)ofatemporal,spectral,andspatialsignature, respec-tively.Unlikematrixdecompositions, thedecompositionof a higher-ordertensorintoasetofsourcesisunique,uptoscalingand permuta-tionambiguities,withoutimposingconstraints(undermildconditions).

ThefMRImatrixissimilarlyapproximatedasasumofrank-1terms. Couplingarisesfromthetemporalsignatures𝐬_𝑟,whichareshared be-tweentheEEGandfMRIfactorization.Afterprocessingthrougha hemo-dynamicsystem(asdescribedinSection2.4.3),eachsource’sBOLD tem-poralsignatureisweightedwithaspatialsignature𝐯_𝐫.

ToaccommodateadditionalstructuredvariationinthefMRIdata, thatis notrelatedtoelectrophysiologicaldynamics, weallowa low-ranktermtothefMRIfactorizationwhichisnotcoupledwiththeEEG factorization.Wehaveempiricallyfoundthatsuchalow-ranktermcan capturestructurednoise,preventingitfrombiasingtheestimationof theparameterswhicharecoupledwiththeEEGfactorization.

ThefullsCMTFmodelisthendescribedas:

=̂+_𝑥 (7) = 𝑅 ∑ 𝑟=1 𝐬𝑟◦ 𝐠𝑟◦ 𝐦𝑟+𝑥 (8) =𝐒,𝐆,𝐌+_𝑥 (9) 𝐘= ̂𝐘+𝐄_𝑦 (10) = 𝑅 ∑ 𝑟=1 𝐾 ∑ 𝑘=1 ( 𝐇𝑘𝐬𝑟)◦(𝐛𝑘∗𝐯𝑟)+ 𝑄 ∑ 𝑞=1 𝐧𝑞◦𝐩𝑞+𝐄𝑦 (11) = 𝐾 ∑ 𝑘=1 ( 𝐇𝑘𝐒) (𝐛T𝑘⊙ 𝐕T ) +𝐍𝐏T+𝐄_𝑦 (12)

=[𝐇1𝐒…𝐇𝐾𝐒]⋅[𝐁T⊙ 𝐕T]+𝐍,𝐏+𝐄𝑦, (13)

where̂ and ̂𝐘arethelow-rankapproximations;𝑥and𝐄𝑦 holdthe

residualsofbothfactorizations;𝐒,𝐆,𝐌describestheCPDmodel com-posedoffactormatrices𝐒∈ℝ𝐼𝑠×𝑅_,𝐆∈ℝ𝐼𝑔×𝑅_,𝐌∈ℝ𝐼𝑚×𝑅_,whichhold

thetemporal,spectralandEEGspatialsignaturesinthecolumns;the HRFmatrices𝐇𝑘areconstructedasin(3)–(6);𝐕∈ℝ𝐼𝑣×𝑅isthefMRI

spatialfactor matrix; 𝐁∈ℝ𝐼𝑣×𝐾 _{is the HRFbasis coefficient matrix;}

𝐍,𝐏isarank-𝑄termtocapturefMRI-onlystructurednuisance;∗ de-notestheHadamardorelementwiseproduct;⊙ denotestheKhatri–Rao product.

Notethatthecoupledpartof𝐘is describedby 𝑅𝐾 nonindepen-dentrank-1terms,orequivalently,by𝐾 rank-𝑅blockterms.Namely, each rank-1 term (𝐇𝑘𝐬𝑟)◦(𝐛𝑘∗𝐯𝑟)describes the convolution of the 𝑟-th source’s temporal signature with the 𝑘-th basis function, after whichaspatialloadingwithvector(𝐛𝑘∗𝐯𝑟)isperformed;inallROIs,

there is one source-nonspecific relative weight for each basis func-tion𝑘(capturedin𝐛𝑘),andsource-specificamplitudes(capturedin𝐯𝑟)

(Calhounetal.,2004).

ItisnotouraimtoestimateHRFvariabilityoversources,butrather, forthesakeofeasierinterpretation,toestimateonlyvariabilityover patientsandROIs.Hence,tolimitthedegreesoffreedom,theHRFin everyROIdoesnotdependon𝑟,butissharedbetweenallsources,as inMaknietal.(2008),Vincentetal.(2010),Pedregosaetal.(2015). Thisis expressedbythetheKhatri–Rao productin(12)–(13), which formsaconstraintthathasearlierbeenusedtorobustifyGLM parame-terestimation(Pedregosaetal.,2015).I.e.,therearenot𝑅𝐾𝐼𝑣spatial

coefficients,but(𝑅+𝐾)𝐼_𝑣,i.e.,𝐾basisfunctionweightsand𝑅source amplitudesin all𝐼𝑣ROIs.Inthisway,theKhatri–Rao structurealso

breaksthecurseofdimensionalityinthefMRIdecompositionifeither thenumberofsources𝑅orthenumberofbasisfunctions𝐾 ishigh(or both).

ThemodelisdepictedinFig.1,omitting𝐍,𝐏tonotoverloadthe diagram.

Weestimateallparametersofthemodelin(8)and(11)byiteratively minimizingthecostfunction𝐽,composedoftwodatafittingtermsand tworegularizationtermsasin(Acaretal.,2014):

𝐽(𝐒,𝐆,𝐌,𝐁,𝐕,𝜽)=𝛽_𝑥||−̂||2_𝐹+𝛽_𝑦||𝐘− ̂𝐘||2_𝐹 +𝛾_𝑥||𝝀_𝑥||₁ + 𝛾_𝑦||𝝀_𝑦||₁ (14) s.t.𝐇𝑘=(𝐡𝑘)=((𝜽(𝑘))) 𝝀𝑥=[𝜆_𝑥,1 … 𝜆𝑥,𝑅 ]T 𝜆𝑥,𝑟=||𝐬𝑟||2⋅ ||𝐠𝑟||2⋅ ||𝐦𝑟||2 𝝀𝑦=[𝜆𝑦,1 … 𝜆𝑦,𝑅]T 𝜆𝑦,𝑟=∑𝐾𝑘=1||𝐛𝑘∗𝐯𝑟||2, (15)

wherethesquaredFrobeniusnorm||||2

𝐹 ofatensoristhesumof

itssquaredelements;||𝐚||2 and||𝐚||1denotetheEuclideanor𝓁2-norm

andthe𝓁1-normorsumoftheelements’absolutevalues ofavector 𝐀,respectively;𝛽𝑥,𝛽𝑦,𝛾𝑥 and𝛾𝑦 arepositive weights;𝝀𝑥and𝝀𝑦 are

vectorswhichholdtheamplitudeswithwhicheachsourceisexpressed intheEEGandfMRIdata,respectively.ThesquaredFrobeniusnorms oftheresidualspromoteagoodfitofthelow-rankapproximationsto thedata,while the𝓁1-regularizationtermspenalizeexcessivesource

amplitudesandpromoteaparsimonious4_{model,similartothe} group-LASSOmethod(Acaretal.,2014;YuanandLin,2006).Atthesametime, thelatterpenaltyalsotendstopreventtheoccurrenceof degenerate terms(Bro,1997).Weminimize(14)usingtheStructuredDataFusion

4 _{Thesparsity-promotingpropertiesoftheLASSOpenaltyaremostusefulin}

thecontextofanunderdeterminedsystem,withmorecoefficientsthan observa-tions,e.g.indictionarylearning.Here,theproblemisheavilyoverdetermined, andwedonotexpectthattheamplitudes𝝀𝑥an𝝀𝑦goexactlytozero.However,

the𝓁1-penaltyisstillausefulheuristictoavoiddegeneratecomponentsinthe

EEG’sCPdecomposition.

frameworkinTensorlab(Sorberetal.,2015;Vervlietetal.,2016),using aquasi-Newtonmethodbasedonalimited-memoryBFGSalgorithm,for 50independentinitializations(seeAppendixAfordetailsregardingthe optimizationprocedureandparameters).Afterconvergence,eachsetof estimatedfactorsneedstobecalibratedtoremovecertainambiguities, andmodelselectionmustbeperformedtopickthebestsolution,with anappropriate𝑅(seeAppendixBfordetails).

2.6. Statisticalinference

Wecreatestatisticalnonparametricmaps(SnPMs)oftheobtained spatialsignatures𝐯𝑟todeterminewhichROIssourcesaresignificantly

(de)activated in relation tothe foundsources (Nichols andHolmes, 2002; Waitesetal.,2005). Namely,underthenullhypothesisof no significantBOLDeffectrelatedtotheEEGdynamics,thefMRIdatamay betemporallyreshuffledwithoutasignificantlossoffittotheEEG dy-namicsin𝐬𝑟.Tothisend,weusepermutation-basedinference,inwhich

thespatialsignatures𝐯𝑟arecomparedagainsttheirempiricallyderived

distributions,whichareobtainedviaresamplingofthefMRIdatawhile freezingtheotherestimatedsCMTFfactors.Toaccountforserial corre-lationsinthefMRItimeseries,weusetherobustwavelet-based resam-plingapproachinBullmoreetal.(2001)toensureexchangeabilityand topreservespatiotemporalcorrelationstructureoftheoriginaldatain theproducedsurrogatedatasets.ForeachfMRIdatasetandeverysCMTF solution,wegenerate𝐿=250surrogatefMRĨ𝐘(𝑙)_{datasetsusingthe}

pro-cedureinBullmoreetal.(2001).Weresampleonlytheadjusteddata

𝐘−𝐍𝐏T,i.e.,afterremovingthecomponentswhich modelvariation specifictothefMRIdata.Weperforminferenceonapseudot-statistic, whichwecomputeforeveryROIandforeverysourceasfollows:

1. constructalocal‘designmatrix’withallestimatedtemporal signa-turesasin(3):𝐃𝑖𝑣= [ 𝐳(1) 𝑖𝑣 …𝐳 (𝑅) 𝑖𝑣 ] ,

2. findthenew‘betas’bysolving𝜷(𝑙) 𝑖𝑣 =𝐃

†

𝑖𝑣̃𝐲 (𝑙) 𝑖𝑣 ∀𝑙,

3. convertthebetastoat-statisticpersourcebydividingthembytheir estimatedstandarddeviation(seeFristonetal.,1994; Polineand Brett,2012).

Throughthisprocedure,weobtain𝐿-pointempiricalnull distribu-tionsforeverysourceandeveryROI.Wesetthesignificancethresholdas tocontrolthefamilywiseerror(FWE)rateat𝛼 =0.05,accordingtothe maximumstatisticprocedureoutlinedinNicholsandHayasaka(2003). Thatis,foreverysource𝑟,weformtheempiricaldistributionofthe max-imalt-statisticoverall𝐼_𝑣ROIs,anddeterminesource-specificthresholds

𝑇(𝑟)

(1−𝛼)asthe95%-percentile(totestforactivation)and𝑇

(𝑟)

(𝛼)asthe

5%-percentile(totestfordeactivation).Finally,weobtainstatisticalmaps forallsources𝑟byapplyingthesethresholdstotheoriginalspatial sig-natures𝐯𝑟,whichcanbeconsideredasthebetasoftheunshuffleddata.

Furthermore,wecreateamapoftheHRFvariabilityoverROIs.For everyROI,weassesshow‘unusual’thelocalHRFis,bymeasuringits calibrateddistanceinHRFspacetoallotherROIs’HRFs.Weusetwo metricstoquantifythis(seeAppendixCfordetailsonthecomputation). 1. Extremityiscomputedasoneminustheaverageoftheabsolute val-uesofthecorrelationsbetweenaHRFwaveformandallotherHRFs’ waveforms.

2. EntropyoftheHRFwaveformiscomputedasthenegativelogarithm oftheconditionalprobabilityoftheHRF.

Bothforthepseudot-mapsasfortheHRFextremityandentropy maps,wefurthermorelimittheinspectiontothe20ROIswiththe high-estvalues,ifapplicable.

Anend-to-endoverviewofourpipeline,fromdatapreprocessingup untilstatisticalinference,isdepictedinFig.2.

2.7. Modelperformance

WeuseseveralmetricstoquantifythequalityoftheobtainedsCMTF solutions.

S. Van Eyndhoven, P. Dupont, S. Tousseyn et al. NeuroImage 228 (2021) 117652

Fig.2. InterictalEEGandfMRIdatacanbeanalyzedviastructuredcoupledmatrix-tensorfactorizations(sCMTF),whichrevealsbothspatiallocalizationofinterictal discharges(spikes),andalsolocalizeddeviationsinneurovascularcouplingbetweenelectricalandBOLDfluctuations.(a)fMRIandEEGdataarefirstseparately preprocessed(yellowblock).ThefMRIdata(toprow)arestructuredasatimepoints×regionsofinterest(ROIs)matrix,afterBOLDtimecoursesareaveragedwithin predefinedordata-drivenparcels.TheEEGdata(bottomrow)arestructuredasachannels×timepoints×frequenciestensor,afterthesignalsareenhancedviaa multi-channelWienerfilter(MWF)whichiscalibratedbasedonspikeannotations,andsubsequentlyundergoatime-frequencytransform.(b)ThesCTMFoftheEEG andfMRIdata(blueblock)revealstemporally,spatiallyandspectrallyresolvedcomponents,andcapturesspatiallyvaryinghemodynamicresponsefunctions(HRFs) (cfr.Fig.1).WeshowtheEEGtemporal,spatialandspectralsignaturesinFigs.4aand6a,andtheHRFsinFigs.4band6b,fortwoselectedpatients.Toinitialize thesCMTFfactors,firstacanonicalpolyadicdecomposition(CPD)oftheEEGtensoriscomputed,fromwhichtheremainingfMRIfactorsareinitialized.Thefull sCMTFmodelisthencomputed𝑁 times,fromthese𝑁 differentinitializations,andthestabilityoftheresultingfactorsoverrunsisassessed.(c)Statisticalimages arecreatedforthepatient’sdataandthecorrespondingsCMTFfactors(greenblock).FromthesCMTFfactors,thespike-relatedcomponentispickedastheone withthehighesttemporalcorrelationtothefilteredEEGsignals’sbroadbandpowerenvelope.Astatisticalnonparametricmap(SnPM)ofthisinterictalspike-related componentiscreated,revealingco-activatedROIsinapseudo-t-map(red).ForeveryROI,theentropy(andalsotheextremity)oftheHRFiscomputedbyassessing itslikelihoodunderthedistributionofallotherROIs’HRFs,andamapofthismetricisconstructed(blue)toreveallocalizedHRFabnormalities.Bothmapscan beusedtoformahypothesisonthelocationoftheepileptogeniczone,asweshowininFigs.5and7forthetwoselectedpatients.Inthispaper,wevalidateour techniqueonasetofpatientsforwhichtheoutcomeisknown.

Wecomparethestatisticalmapswithagroundtruthdelineationof theictalonsetzone(IOZ)toassesstheirconcordance.Thisgroundtruth isthemanuallydelineatedresectionzoneforpatientsthathad under-gonesurgicaltreatmentandthatwereseizure-freeafterwards(Anetal., 2013;Grouilleretal.,2011; vanHoudtetal.,2013; Thorntonetal., 2010;Zijlmans etal., 2007),orotherwisethehypotheticalresection zone,basedonconcordantevidencefrommultiplemodalitiesotherthan EEG–fMRI(cfr.Section2.1),forpatientsthatwereineligibleforor re-fusedsurgery(Tousseynetal.,2014a).Thesensitivityfordetectingthe IOZisthencomputedasthefractionof‘truepositive’cases,whichare determinedbythepresenceorabsenceofsignificantactivationclusters whichoverlaptheIOZinthespatialsignatures𝐯𝑟.Followingthe

reason-inginTousseynetal.(2014a),wedonotconsidersignificantlyactive voxelsorregionsoutsideofthedelineatedIOZasfalsepositives. Ac-knowledgingepilepsyasanetworkdisorder,suchactiveregionsmight reflectseizureorIEDpropagation,despitenotbeinginvolvedintheir generation.

Furthermore,wehypothesizethatthespatialvariationoftheHRF overthebrainmightrevealadditionallocalizinginformation regard-ingtheIOZ,i.e.,basedon considerationsexplainedin Section1, we

assumethattheHRFinorneartheIOZmightbedistortedcomparedto nonepilepticbrainregions.Wetestthishypothesisbyassessingwhether thoseregionscorrespondtohighvaluesintheHRFentropyandHRF ex-tremitymaps(cfr.Section2.6).

Additionally,weinspectthespectral,spatialandtemporalEEG sig-naturesoftheextractedsources,andwemeasurewhetherthespatial fMRIsignaturesbearanysimilaritytoknownnetworksofresting-state humanbrainactivity(Shireretal.,2012).

3. Experiments

3.1. Patient-specificmodelselection

TableB.3compilestheresultsofthemodelselectiondescribedin

AppendixB.Foreachpatient,weselectthesetofsCMTFfactorsofrank

̂𝑅,whichbestfulfillthecriteria.Inallcases,wefoundatleastonesucha solution,includinganIED-relatedcomponentwithinthatsolution.Note thatsometimesmodelswithdifferent𝑅mightscorewellondifferent (subsetsofthe)criteria,sotheselectionoftherankisinevitably am-biguous.Inthenextsection,weanalyzetheindividualsetofresultsfor

Fig.3.Goodnessoffitofeachpatient’sEEGtensor andfMRImatrix𝐘,forvaryingchoicesoftherank𝑅

in thesCMTF. Naturally,theEEG approximation er-rordecreasesmonotonicallyforincreasingrank (intra-patient).ForthefMRIdata,thefitalreadyplateausfor verylow𝑅.Thisisduetothepresenceofadditional, uncoupledcomponents𝐧𝑞◦𝐩𝑞inthefMRIfactorization,

whichcanabsorbsomeofthevariancewhenthe num-berofcoupledcomponentsislow,butwhichlosetheir relevanceathigherranks.

eachpatient,basedontheselectedrank,andweanalyzethesensitivity oftheresultstothechoiceof𝑅.

Weshowthegoodnessof fitoftheestimatedfactorsfortheEEG tensor andthefMRImatrix𝐘inFig.3.Duetothenormalization stepswhichhavebeenappliedtothedata(cfr.Section2.2),thesCMTF operatesinaregimeofmoderatelyhighrelativeapproximationerrors.

3.2. Spatio-temporo-spectralprofilesofinterictaldischarges

Weanalyzeforeachpatientthesourceswhichhavebeenestimated viathesCMTFmodel.Wediscusstheresultsoftwopatientsindetailin thenextsubsections,andincludecompleteresultsforallotherpatients inthesupplementarymaterial.

Everytime,weshow(1)thethresholdedpseudot-mapsofthe IED-relatedsourceinthefMRIdomain,bothforsignificantactivationasfor significantdeactivation;(2)mapshighlightingtheROIsof highHRF entropyandextremity;(3)thetemporalprofile(time-varyingpower)

𝐬𝑟,spatialprofile(topography)𝐦𝑟andspectralprofile𝐠𝑟ofeachsource

intheEEGdomain;(4)theHRFwaveformsinthedifferentROIs,and theHRFbasisfunctionsatconvergenceofthealgorithm.

Weplotmaximally800softhetemporalsignatures,toensure read-ability.Foreaseofcomparison,wealwaysoverlaythebroadbandMWF envelope(withanarbitraryverticaloffsetforvisualizationonly),which isthereferencetimecourse𝐬refforselectingtheIED-relatedcomponent

(cfr.B.3).Forconsiderationsofspace,wegenerallyonlyshowthemaps ofthefMRIspatialsignature𝐯𝑟fortheIED-relatedcomponents,but

dis-cussthemapsofothercomponentswhenrelevant.Weshowfiveaxial slicesofeachmap:ineachcase,weshowtwoslicesnearthehighest andlowestvoxelsoftheIOZorsignificantregionsofthefMRIspatial signature(whicheverliesfurthest);ifapplicable,themiddlesliceisthe cross-sectionwithmostoverlapbetweenIOZandspatialsignature,and thetworemainingslicesliehalfwaybetweenthissliceandtheextremal slices;otherwiseallthreebulkslicesarechosenwithequalspacing be-tweentheextremalslices. Wecross-validate themapsagainstknown restingstatenetworks(RSNs)ofhumanbrainactivityfromtheStanford atlas(Shireretal.,2012).

Westressagainatthispointthatasubsetoftheresultsisproneto er-rorsduetoimperfectsignnormalization(cfr.B.1).Whileitisrelatively straightforwardtounambiguouslydeterminethe‘right’signoftheEEG signatures,thisismorechallengingforfMRI. Thatis,frequently,the polarityoftheHRFwaveformis ambiguous,andmakingthe‘wrong’ choiceinavoxel𝑖𝑣(i.e.,theHRFhastheoppositeeffectofthetrue

physicalcerebralbloodflowchange)immediatelyleadstothewrong signofthespatialcoefficientsin𝐯_𝑟intherespectivevoxel,andtheir pseudot-values,forallsources𝑟.Totracktheoccurrenceofthis

fore-seenfailuremode,wealsoinvestigatethesignificantdeactivationsof thesources5_.

NotethatwedesignedtheHRFvariabilitymetricssothattheyare im-munetothepolarityoftheHRFs.Hence,anyhighscoreoftheHRF met-ricscanbereliablyinterpreted.Foreachcase,weseparatethetwenty waveformswiththehighestentropyscores,andreporthowmanyof thosearefoundinROIsthatoverlapwiththeIOZ,alongwiththe prob-ability(intheformofap-valuefromabinomialdistribution)thatthis wouldoccurbyrandomlysamplingasmanyROIs(underagiven frac-tionofbrainthatiscoveredbytheIOZ).Hence,thismetricisanalogous tooneminusthefalsediscoveryrate(FDR).

3.2.1. Patient3

Weanalyzethesolutionwith ̂𝑅=2sources,andshowtheresultsin

Figs.4and5.BesidesoneclearIED-relatedsource,thereisoneother sourcethatissubstantiallycorrelatedtothereferencetimecourse,but withahomogeneousdistributionover theheadandanunclear spec-trum. Thismaysignify thattheIEDs donot follow exactlyarank-1 structure inthespectrogram,andthattheymaybe nonstationaryin time orspace (cfr.theargumentmade fornonstationary seizuresin

Hunyadi etal.,2014). Thesecondsource’spseudo t-maphad signif-icantlyactiveareas symmetricallyintheleftandrightparietallobe, muchmorefocalizedthantheEEGtopography.IntheEEGtimecourses, wefoundindeedIED-likewaveformsatthetimesofthepeaksinthe temporalsignature.Hence,wesuspectthatbothsourcesmayreflectthe onsetandpropagationoftheIEDstootherareas,respectively.Fiveout ofthetwentyROIswithhigh-entropyHRFsoverlappedwiththeIOZ, andasignificantfindingisthatseveralofthemarehighlynoncausal, i.e.,withapositivepeakbeforezeroseconds.Fig.5confirmsthis,and alsoshowsthattheIED-relatedsourceissignificantlyactiveindifferent ROIsoftheIOZ.

3.2.2. Patient10

Weanalyzethesolutionwith ̂𝑅=5sources,andshowtheresultsin

Figs.6and7.ThereisaclearIED-relatedsource,andalsoan artifac-tualsourceat±33Hz,whichisalsopresentinotherpatients.Duetoits relativelyconsistentoccurrence,wehypothesizethatthisartifactisdue totheMRacquisition.Forexample,itmaybearemnantofagradient

5_{Alternatively, itispossibletousea pseudoF-statistic,e.g. thesquared}

pseudot-value,tobypassthesigncorrectionaltogether.Thedownsideofsuch anapproachisthatitisthenimpossibletodistinguishactivationand deactiva-tion,whichmaybemeaningful.

S. Van Eyndhoven, P. Dupont, S. Tousseyn et al. NeuroImage 228 (2021) 117652

Fig.4. (a)Intheselectedsolutionforpatient3(̂𝑅=2),bothsourceshaveatemporalsignaturethatcorrelatedstronglytothereferenceIEDtimecourse.Thefirst sourcemodeledthemainonsetofIEDsandwaslow-frequencyandtopographicallyfocal,whilethesecondsourcewasspatiallyandspectrallydiffuseandcaptured thepropagationofIEDstoremoteareas(cfr.Fig.8).(b)FiveoutofthetwentymostdeviantHRFswerefoundinsidetheictalonsetzone(boldlines,𝑝<10−4_).These

HRFshadmainpeaksbefore0s,i.e.,theyledtoBOLDchangesbeforethecorrespondingEEGcorrelateoftheIEDwasseen.

artifactwhichisnotadequatelyremovedfromthedataofsome chan-nels,cfr.theobservationmadeinMarečeketal.(2016).Surprisingly, thissourceissignificantlyactiveinanextendedareaintheoccipital lobe,overlappingwiththevisualnetwork.BothHRFmetricsreached extremevaluesatsome (distinct)ROIswithin theIOZ.Thepseudo-t mapoftheIED-relatedsourceshowssignificantlyactiveROIsthatare concordantwiththeIOZ,anddeactivationofalargepartofthedefault modenetwork.Furthermore,theIED-relatedsource’sEEGtopography isveryconsistentwiththeclinicaldiagnosis.Thefourthsourceisactive inthedefaultmodenetwork,predominantlyinthe𝛼 band(cfr.Fig.9).

Thefifth sourcehadanunclear spectrum,butits temporalsignature correspondstotheoccurrenceofhigh-amplitudeIEDs.Itspseudot-map showswidespreadactivationsoverthebrain,whichdidnotincludethe IOZ.WeexpectthatthiscomponentcapturesthepropagationofIEDs, afteronsetneartheIOZ,similarlytopatient3.

3.2.3. Summaryofallpatient’sresults

Weprovideanoverviewoftheresultsw.r.t.IOZdetectioninTable2. Allresultstakentogether,thesCMTFresultsallowacorrectdetection of theIOZbasedonthesignificantIEDactivation(10/12cases)and

Fig.5.ThestatisticalnonparametricmapsoftheIED-relatedcomponent(toptworows)andHRFentropy/extremitymaps(bottomtworows)ofpatient3show concordancewiththeictalonsetzone(IOZ).EspeciallytheregionsofsignificantIEDactivationwereaccurate,butalsofiveoutofthetwentyregionswiththemost deviant(highestentropy)HRFswerefoundintheIOZ(cfr.Fig.4b).Thegroundtruthictalonsetzoneishighlightedindarkgraywithawhitecontour.ROIswith highvaluesforbothHRFvariabilitymetricsarecoloredinorange.

Table2

ThesCMTFleadstothreetypesofspatialinformation,whichcanbecross-validatedagainstthegroundtruthIOZ,asdefinedinSection2.7andsummarizedforall patientsinTable1:(1)theEEGtopography𝐦IEDoftheIED-relatedcomponent;(2)thesignificantlyactivatedanddeactivatedROIsinthefMRIspatialsignature

𝐯IED;(3)theROIswithstronglydeviatingHRFwaveforms,asmeasuredviaentropyandextremity.SincetheEEGtopographyhasaverylowspatialresolution,

anddependsontheattenuationpropertiesofthetissueaswellastheorientationoftheneuralsourcesinthecortex,weonlyexpectpartialsimilaritytotheIOZ’s spatialfocus;hence,weusetheterm‘consistent’ratherthan‘concordant’.Thepatientswhohadagoodoutcomeaftersurgery(patients5,7and11)hadahigher concordancebetweenthethreetypesofspatialcluesthanpatientswithapooreroutcome(patients1,2and6).

patient selected solution

EEG topography consistent with IOZ?

spatial signature 𝐯 IED concordant with IOZ?

HRF variability metrics complementary to 𝐯 IED ?

20 highest-entropy ROIs

ID ̂𝑅 activation deactivation entropy extremity # in IOZ (p-value)

p01 6 ✓ ✓ 1 (0.34) p02 3 ✓ ✓ ✓ 1 (0.59) p03 2 ✓ ✓ ✓ 5 ( < 10 −4 ) p04 4 ✓ ✓ ✓ ✓ 2 (0.32) p05 5 ✓ ✓ ✓ ✓ ✓ 6 ( < 10 −3 ) p06 2 ✓ 0 / p07 4 ✓ ✓ 1 (0.57) p08 2 ✓ ✓ 0 / p09 2 ✓ ✓ ✓ ✓ 0 / p10 5 ✓ ✓ ✓ ✓ 3 (0.02) p11 2 ✓ ✓ ✓ ✓ ✓ 4 (0.01) p12 2 ✓ ✓ ✓ 7 ( < 10 −3 )

significantIEDdeactivation(6/12cases).Formanypatients,someof theROIswiththehighestHRFentropy(9/12cases,ofwhich8were complementarytotheSnPM)andhighestHRFextremity(8/12cases) alsooverlappedwiththeIOZ,whichwasshowntobe(very)unlikely duetochance.Allcasesarecoveredbyatleastoneofthemetrics,and allpatientsbesidespatient6hadatleasttwometricsprovidingcorrect andcomplementarylocalizinginfoontheIOZ.Fornearlyallcases,the IED-relatedcomponent’stimecoursewashighlycorrelatedtoa refer-enceIEDtimecourse,anditsspectrumwasplausible.Inmany,butnot allcases,thiscomponent’sEEGtopographywasalsoconsistentwiththe locationoftheIOZ,thoughthisnotionisslightlyfuzzybecauseofthe verydifferentspatialdomainsofEEGand(f)MRI—hencewedonotuse theterm‘concordant’.Analysisof thespatial,spectral,andtemporal signatures,incombinationwithinspectionofthefilteredEEGsignals, revealstheidentityof RSNoscillationsand/orartifactsin the major-ityofcases.Forseveralpatients,wefoundsourcesthatareactivein anarrowspectralbandnear33Hz.Whiletheselikelyreflecta techni-calartifactastheresultoftheMRacquisition,wefoundno concomi-tantchangesatthisfrequencyintheEEG.Thismaybetheresultofthe normalizationprocedurewhichweappliedpriortothedecomposition: sinceeveryfrequencybinwasgivenequalimportance,even

unnotice-ablebutstructuredfluctuationsathigherfrequenciesmaybecaptured inacomponent.

3.2.4. Sensitivitytomodelselection

Formanypatients,selecting ̂𝑅isambiguous,sincemorethanone solution (with different 𝑅) score well on some of the criteria (cfr.

TableB.3).Therefore,weanalyzetheimpactofthechoiceof𝑅onthe sCMTFresults.Foreachpatient,weselectthesolutionwiththerank which isnextinline,i.e.,whichwouldbe asecondbest(orequally good)choice,basedonthesamecriteria.Thisisthesolutionwith𝑅=1 forpatient12,𝑅=2forpatients1,2,5and7,𝑅=3forpatients3,4, 6,8,9and10,and𝑅=4forpatient11.Forpatients1,6and8,the resultsdeterioratedrastically,asnometriccorrectlylocalizestheIOZ. Forpatient11,noROIwithintheIOZissignificantlyactivateddueto IEDsanymore,buttheHRFmetricsarestillinformative.Theresultsfor patients9and12improve,sinceallmetricsarenowsensitivetotheIOZ. Fortheotherpatients,thesituationstaysmoreorlessthesame,i.e.,the samemetricsarevaluableforIOZlocalization.However,themaximum valueunderthedifferentmetricsisgenerallyattainedatdifferentROIs comparedtotheinitiallyselectedmodel.

S. Van Eyndhoven, P. Dupont, S. Tousseyn et al. NeuroImage 228 (2021) 117652

Fig.6. (a)ThesCMTFsolutionwith ̂𝑅=5sourceswasselectedforpatient10.Onesource’stemporalsignatureishighlycorrelatedwiththereferenceIEDtime courseandisidentifiedastheIED-relatedsource,whichhasacharacteristiclow-frequencybehaviourandwithafrontotemporaltopography,consistentwiththe IOZlocation.Thesecondsource,whichhasverynarrow-bandpoweraround±33Hz,likelycapturedanartifactoftheMRacquisition.Thefourthsourcecaptured𝛼

activityinthedefaultmodenetwork(cfr.alsoFig.9).(b)ThreeoutofthetwentymostdeviantHRFswerefoundinsidetheictalonsetzone(boldlines,𝑝=0.02).

4. Discussion

AnovelEEG–fMRIdatafusionframework

We have proposed an integrated and structured coupled matrix-tensorfactorization(sCMTF) framework,whichcanbe usedtomake inferencesonthelocalizationoftheictalonsetzoneinrefractoryfocal epilepsybasedonsimultaneousEEGandfMRIrecordings.Ourapproach aimstoperformblindsourceseparationoftheneuralactivityrelatedto interictalepilepticdischarges(IEDs),andtocharacterizeitinthe spa-tial,temporal,andspectraldomain.Tothisend,wedevelopedapipeline consistingof(1)semi-automatedEEGenhancementbasedon

annota-tionsoftheIEDs;(2)modality-specificpreprocessingandtensorization steps,whichleadtoathird-orderEEGspectrogramtensorvaryingover electrodes,timepoints,andfrequencies,andanfMRImatrixwithBOLD timecoursesforapredefinedsetofregionsofinterestorparcels;(3) coupledmatrix-tensorfactorizationoftheEEGtensorandfMRImatrix alongthesharedtemporalmode,whileaccountingforvariationsinthe localneurovascularcoupling;(4)automatedselectionofarobust,and relevantIED-relatedcomponent,andnonparametrictestingtoinferits spatialdistributioninthebrain.

Wehavestressedtheimportanceofandaccountedforthevariability of thehemodynamicresponsefunction (HRF)over differentpatients

Fig.7.ThestatisticalnonparametricmapsoftheIED-relatedcomponentandHRFentropy/extremitymapsofpatient10showconcordancewiththeictalonsetzone (IOZ).IEDoccurrenceswereassociatedwithsignificantlyactive(red)regionsinandneartheIOZ,andatthesametimewithadeactivation(blue)inapartofthe defaultmodenetwork.Threeoutofthetwentyregionswiththemostdeviant(highestentropy)HRFswerefoundintheIOZ(cfr.Fig.6b).Thegroundtruthictal onsetzoneishighlightedindarkgraywithawhitecontour.ROIswithhighvaluesforbothHRFvariabilitymetricsarecoloredinorange.

Fig.8.Thesecondcomponentinpatient3likelycapturedthepropagationofIEDsfromtheirritativezone,givenitsrelativelylargecorrelationtotheMWFenvelope (cfr.Fig.4a).Thegroundtruthictalonsetzoneishighlightedindarkgraywithawhitecontour.

Fig.9. Thefourthcomponentinpatient10seemedtopickupactivityintheDefaultModeNetwork(DMN),predominantlyinthe𝛼 band(cfr.Fig.6a).

andbrainregions,byequippingtheCMTFwiththerequiredexpressive powerviaasetofadaptivebasisfunctions.Moreover,afterestimating theEEGandfMRIfactorsignatures,aswellastheHRFparameters,we havecomputeddifferentsummarymetrics(entropyandextremity)that measurethelocaldevianceofaROI’sHRFcomparedtootherHRFsin thesamebrain,andhavecross-validatedthespatialmapofthesemetrics againstthegroundtruthlocalizationoftheictalonsetzone.

ThesCMTFpipelineprovedtobesensitiveindetectingtheIOZinall twelvepatientsinthisstudy.Thestatisticalnonparametricmap(SnPM) ofthespatialsignatureoftheIED-relatedcomponent’sactivation, ob-tainedwith thesCMTF,is thebestbiomarker, whichis in linewith thetraditionalEEG-correlatedfMRIapproach(Lemieuxetal.,2001). Inthelargemajorityof patients,severalof thesesignificantlyactive ROIsoverlappedwiththeIOZ.Whenusedinconjunctionwiththe IED-relatedSnPM,alsotheHRFentropyandextremity,asmeasuresofhow unlikelyanHRFiswithinaspecificsetofotherHRFs,arepromising biomarkers,whichidentifiedregionsoftheIOZthatwere complemen-tarytothosealreadyfoundbytrackingsignificantIEDactivationfor

nineoutoftwelvepatientsinthisstudy.Inroughlyhalfofallcases, wealsofoundregionswithintheIOZthatsignificantlydeactivatedin associationtoIEDs.Thepatientswhohadagoodoutcomeaftersurgery (patients5,7and11)hadahigherconcordancebetweenthethreetypes ofspatialclues(EEGtopography,fMRI(de)activation,HRFvariability) thanpatientswithapooreroutcome(patients1,2and6).Whilethe numberofpatientsistoolowtodrawconclusions,thisobservation sup-portsthehypothesisthatthedegreeofsuchconcordancemighthelpto predictsurgicaloutcome.Initscurrentform,thepipelinestillpredicted toomanyROIsthatdidnotoverlapwiththeIOZ.Thismightinpartbe duetoIEDpropagationthroughouttheepilepticnetwork,aspostulated inTousseynetal.(2014a),butislikelyalsoaresultofinherent limita-tionsofthemodel.Hence,theoutputofthisanalysisisonlytobeused inconjunctionwithothermodalities(e.g.SISCOM)duringpresurgical evaluation,inorder toassesscross-modalityagreement,asisalready commonforEEG–fMRI.

Wepreviouslystudiedthesamepatientcohortusingclassical EEG-correlatedfMRIanalysis,usingdifferenttypesof EEG-derived

regres-S. Van Eyndhoven, P. Dupont, S. Tousseyn et al. NeuroImage 228 (2021) 117652

sors,inVanEyndhovenetal.(2019a). There,weconcludedthatthe time-varyingpowerofMWF-enhancedEEGwaspreferableover regres-sorsbasedonICA,EEGsynchronization,orsimplestickfunctions, al-lowingasensitivityoftheIEDactivationmapforIOZlocalizationof elevenoutoftwelvecases.Comparedtothispriorstudy,thecurrent pipelineyieldedanadditionalcorrectdetectionforpatient11,butfailed tocorrectlypredictROIsthatoverlapwiththeIOZforpatients1and 12.However,forpatient12,theHRFentropyandextremitymetrics wereveryinformative,inthattheypredictedmanyIOZROIscorrectly. Hence,theadditionalflexibilityofthis pipelinewas probablykeyto theIOZdetectionforpatient11.Yet,thetwo‘misses’fortheothertwo patientsindicatethattherobustnessofthecurrentdesignstillneeds im-provement.Forcomparison,weincludeinthesupplementarymaterial thestatisticalmapswhichweobtainedviatheanalysisinVan Eynd-hovenetal.(2019a).Itisworthnotingthatbothpatient1andpatient 11hadonly fewIEDs(15and6,respectively),whichmakes mining oftheir EEG–fMRIdatafor BOLDactivationdifficultapriori( Salek-Haddadietal.,2006;2003;Zijlmansetal.,2007).Hence,alsoforour moreflexiblemethod,theprobableyieldscaleswiththenumberofIEDs duringrecording.Weinspectedthe20HRFsandROIswiththe high-estextremityandentropy.Hence,itisinevitablethatsomeormostof theseROIsarenotwithintheIOZ,ortheIOZmightevennothavea deviantHRF.StandaloneHRFmetricswouldhencehaveahighfalse discoveryrate,eventhoughforseveralpatients,thehighproportionof IOZ-coveringROIsamongthe20selectedROIswasveryunlikelydueto chance(asmeasuredwithp-values).Still,itisnotalwaysthecasethat theHRFinsidetheIOZismeasurablydifferentthanintherestofthe brain,inwhichcasethesemetricswouldnotbeverysensitive.

However,theROIsthatwerehighlightedbytheHRFmetricswere oftendistinctfromtheROIsidentifiedassignificantlyactivatedtothe IEDs.Hence,theSnPMsoftheIEDsandtheentropymetricsprovidevery complementaryinformation,andwhenanalyzedjointly,theymayinfer thelocationof theIOZwithmorecertainty,byselecting brainareas wherebothIED-relatedandHRF-relatedmetricshaveahighvalue.

Weenvisionthatthisapproach,withminormodifications,mayalso beusedtoanalyzeresting-stateEEG–fMRIactivity,sincealsothe estima-tionofextendednetworksofcorrelatedspontaneousBOLDfluctuations canbebiasedbyspatiallyvaryingHRFs.Sinceinsuchacase,noIEDs arepresent,EEG-enhancementlikewehavedoneinthisstudywouldno longertakeplace.However,anMWFmaystillbeusedtocleanupthe EEG,e.g.byannotatingartifactualperiods,whichcanberemovedfrom thedatabytheMWFinitsdualform(oranothertool)(Somersetal., 2018).SimilartothemethodinMartínez-Montesetal.(2004),the cur-rentmethodcouldthenallowtoextract RSNswhicharereflectedin bothEEGandfMRIdata,giventhatsufficientlymany components𝑅

areextracted.

HRFsvarystronglyoversubjectsandbrainregions

Thereweresubstantialdifferencesin(estimated)neurovascular cou-plingoverpatientsandbrainregions,asexpected.Sinceweused ‘regu-larized’basisfunctions,whichareparametrizedassmoothgamma den-sityfunctions,theresultingHRFsgenerallyhadaplausibleshape. How-ever,in somecaseswefoundnonsensicalshapes, inwhich,e.g.,the waveformahadthesamepolarityoverthewholetimecourse, poten-tiallywithabimodal shape(cfr.patient4).This servesasahumble remindertonotblindlytrusttheoutputtedHRFs(orotherfactor sig-natures,forthatmatter).Whilewehaveempiricallyverifiedthatthe optimizationalgorithmconvergesproperlytothetruefactorsignatures andHRFsforsyntheticdataundermildconditions,thereisno guaran-teethatthisholdstrueforreal-lifedata,whichareordersofmagnitude morecomplex,sothatalineargenerativemodellikethesCMTFmay notbesufficienttodescribetheinterplaybetweenEEGandfMRI. More-over,theproperbehaviorofthesCMTFestimationdependsoncareful preprocessing,andonaproperselectionofhyperparameters(incasu:a goodvalueforthenumberofsources ̂𝑅).Hence,manualinspectionof thedataqualityandthesolutionisstillrequired.Eveniftheestimated

HRFsorfactorsignaturesmaynotfullyreflectthe‘correct’underlying physicalphenomena,wehavedemonstratedthattheyofferactionable information.Notintheleast,viasummarizingmetricssuchasHRF en-tropy andextremity,ouralgorithmmanagestobe reasonablyrobust tosubtlechangesinthewaveform—whichislessofinterestherethan spatialcuestowardstheIOZ.Wereckonthatotherparametrizationsfor HRFthantheonewehaveused,mightbebettersuitedforthetask. These basisfunctionscould evenbe pickedapriori,e.g., froma set ofsensiblegeneratingparameters(Woolrichetal.,2004).Thiswould evensimplifytheoptimizationproblem,sincetheparametervectors𝜽𝑘

nolongerneedtobeestimated,attheexpenseofbeinglessdata-driven. Thealgorithmusedits modelingfreedomtofit‘noncausal’HRFs, whichareaheadoftheEEGbyasmuchas10s.Generally,weindeed foundthatmanyoftheestimatedHRFshadsignificantpositiveor neg-ativeamplitudesalreadybeforetheneuralcorrelatevisibleontheEEG. This isinlinewithrecurrentfindingsonBOLDchangesthatprecede theIEDswhichwereobservedintheEEG(Hawcoetal.,2007;Jacobs etal.,2009;Moelleretal.,2008;Pittauetal.,2011).Westressthatthis noncausalitymayonlybeintheobservation,andnotintheunderlying physicalchainofevents:here,itstrictlymeansthatweobservedBOLD changesinthefMRIdatathatoccurbeforethecorrespondingobserved

neuralcorrelateontheEEG.DespitethefactthatmanyoftheHRFs dif-feredsubstantiallyfromthecanonicalHRF,whichiscausalandpeaks approximately 6safter itsneuralinput,we obtainedgoodresults as wellwiththelatterHRFasanonadaptivemodelforneurovascular cou-pling(VanEyndhovenetal.,2019a).Thesameconclusionwasreached by Caballero-Gaudesetal.(2013)for thecomparisonof the canoni-calHRFandaninformation-theoreticapproachforBOLDmapping.The reasonforthisagreementbetweenthesedifferentmodels—whichdiffer substantiallyintermsofflexibility—islikelythatthecanonicalHRFis positivelycorrelatedtothetrueHRFswhicharefoundinsidetheIOZ, andassuchtheresultingactivationmapsmaystillbesufficiently infor-mative.InourdataandsCMTFresultsthisisindeedthecaseformany patients.

PriorEEGsignalenhancementaidsanalysis

Importantly,ourpipelineheavilyreliesonapriorenhancementof theinterictalspikesintheEEGdata,whichwouldotherwisehaveatoo lowSNRforthesCMTFalgorithmtopickupIED-relatedsources.We employmulti-channelWienerfilters,whichsolelyrelyonthe annota-tionofasufficientamountofIEDsinthedataitself,orinrelateddata (e.g.,datafromthesamepatient,recordedoutsidetheMRscanner). WhilethistaskstillfrequentlyreliesontheskillofhumanEEGreaders andneurologists,advancedautomatedsolutionsforinterictalspike de-tectionareavailable(Scheueretal.,2017;Wilsonetal.,1999).Within eachsolutionofaspecificrank,wepickedtheIEDcomponentasthe onewiththehighestcorrelationwithareferencetimecoursedirectly derivedfromtheenhancedEEG.Someof thepresented resultsmake clearthatthisreferencetimecourseisnotcompletelyfreefrom arti-facts,hencecautioniswarrantedwhenmanyhigh-amplitudeartifacts arestillpresentinthereference.Inthisstudy,however,wehavenot encounteredanyissuesthatseemedtobethedirectresultsofanoisy referenceduringIEDcomponentselection.ForthefMRIdata,wehave carriedoutarelativestrict,butunsupervised‘enhancement’,by regress-ingoutmultiplepotentialconfounds.Hence,itwouldbeworthwhileto performthefMRIcleanupaccordingtoatask-basedorsupervised crite-rionaswell,e.g.,usingICAcombinedwithnoisecomponent identifica-tion(Salimi-Khorshidietal.,2014).

Theinterpretationofcomponents

Overall,thesCMTFpipelinesucceededinextractingmeaningful IED-relatedcomponents,alongsidecomponentsthatmodeled resting-state neuralfluctuationsandphysiologicalandtechnicalartifacts.Thefact thatthesCMTFcanestimatesignaturesandstatisticalmapsformultiple componentsisapowerfuladvantageoverclassicalEEG-correlated anal-ysis.Aswedemonstratedintheexperiments,artifactualinfluencesmay