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

(1)

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

a,∗

_,

_Patrick

_Dupont

b,c

_,

_Simon

_Tousseyn

d

_,

_Nico

_Vervliet

a

_,

_Wim

_Van

Paesschen

e,f

_,

_Sabine

_Van

_Huﬀel

a

_,

_Borbála

_Hunyadi

g

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 aneurologicaldisordersuﬀeredby30%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)signalsforsigniﬁcant

cor-∗_{Corresponding}_author.

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 reduceitssensitivityorspeciﬁcity(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

(2)

Firstofall,theadoptionofarelevantrepresentationofIED occur-rencestoconstructaregressorforthedesignmatrixhasprovenvitalto thesensitivity.ThisaspecthasbeeninvestigatedinRosaetal.(2010), Murtaetal.(2015),Abreuetal.(2018),VanEyndhovenetal.(2019a). Inpreviouswork(VanEyndhovenetal.,2019a),weaddressedthis is-suebypre-enhancingtheEEGsignalsusingaspatiotemporalﬁlterthatis tunedtomaximizethesignal-to-noiseratio(SNR)ofIEDswithrespect tothebackgroundEEG.Wehaveshownthattakingthetime-varying poweroftheﬁlteredEEGleadstoarobustregressor,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 opportunitytomapnotonlyregionswithstatisticallysigniﬁcantBOLD changesinresponsetoIEDs,butalsothespatialmodulationoftheHRF waveformitself, inorder todiscoverregions whereanaﬀectedHRF 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

(3)

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-tialweightingcoeﬃcientsareestimatedalongwithspatial,spectraland temporalsignaturesofcomponents.BypreprocessingtheEEGusingthe data-drivenﬁltersfromVanEyndhovenetal.(2019a),weaimto max-imizethesensitivityinmappingtheinterictaldischarges.Weanalyze whethertheestimatedspatialmodulationoftheHRFwaveformisa vi-ablebiomarkerwhenlocalizingtheictalonsetzone,besidestheBOLD spatialsignaturesthemselves.

2. Methodsandmaterials 2.1. Patient group

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 classiﬁcation1-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. Data acquisition and preprocessing

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-passﬁlteredoﬄinebetween1-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

(4)

atlas,whichconsistsof246parcelsinthegreymatter(Fanetal.,2016). ForeveryROI,oneBOLDtimeseriesisconstructedastheaverageof thetimeseriesofallvoxelswithintheROI.Ifmultipleacquisitionruns (withinthesamerecordingsession)hadbeendone,theEEGandfMRI dataof the diﬀerent runs aretemporally concatenated. Further cus-tomizedpreprocessingstepsaretreatedinSections2.3and2.4. 2.3. Multi-channel Wiener ﬁltering for spatio-temporal EEG enhancement

Inpreviouswork(VanEyndhovenetal.,2019a),wehaveshownthat pre-enhancingtheEEGsignalsusingadata-driven,spatiotemporalﬁlter 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-channelWienerﬁlters(MWF)haserror-correctingcapabilities andproducesanIEDrepresentationthatimprovesthelocalization sen-sitivityofEEG-correlatedfMRI(VanEyndhovenetal.,2019a).

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

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

andsubsequentlycomputingtheﬁltercoeﬃcientsas ̂

𝐖=𝐑−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_̃𝐱_. _Due_to_the_extension_with_lagged_copies_of thesignals,channel-specificfiniteimpulse responsefiltersarefound. Hence,𝐖̂T_̃𝐱_is_a_set_of_{spatiotemporally}_filtered_output_signals,_in_which IED-like waveformsarepreserved whileother waveforms,whichare notspecifictoepilepsy,aresupressed1_.

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

2.4. Higher-order data representation

Topreservetheintrinsicmultiwaynatureofthedata,werepresent thepreprocessedEEG andfMRI asatensorandmatrixrespectively, whicharesubsequentlyfactorizedjointly.Thisapproachdiﬀers from themass-univariatetreatmentinthetraditionalGLM,whereeachvoxel istreatedindividually,andonly‘ﬂattened’EEGtimecoursescanbe en-teredasregressors.Sinceepilepsyismanifestedwithconsiderable vari-abilitybetweenpatients,wehandlethemultimodaldataofeachpatient separately.

1 _To_retrieve_ﬁltered_versions_of_the_original_set_of_channels_only_(and_not_of theirtime-delayembeddedcopies,whichwouldberedundant),onlythemiddle

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

2 _We_observed_in_Van_Eyndhoven_et_al._(2019a)_that_for_most_patients,_the outputSNRsaturatesaroundthisvalue,correspondingtoanintervalof-16ms to+16ms.

2.4.1. Spatio-temporal-spectral tensor representation of EEG

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∈ℝ𝐼𝑠 ×𝐼_𝑔×𝐼_𝑚_is_further_normalized_as_described_in_A.1_,_to equal-izetheinﬂuenceofeachchannelandeachfrequency,andtofocuson relativesignalincreasesordecreases(Marečeketal.,2017;2016) 2.4.2. Spatio-temporal matrix representation of fMRI

TheaverageBOLDtimeseriesarestackedinatimepoints×ROIs matrix𝐘∈ℝ𝐼𝑠 ×𝐼𝑣 _,where𝐼𝑣=246ROIs.WenormalizeeachROI’stime seriesbysubtractingitsmeananddividingbyitsstandarddeviation. 2.4.3. Neurovascular coupling in the temporal mode

EEGandfMRIdataareacquiredsimultaneouslypersubject,andare thusnaturallycoregisteredalongthe‘time’mode.Thisiscapturedina temporalfactormatrixthatiscommonbetweentheEEGfactorization andthefMRIfactorization.However,theelectrophysiologicalchanges thatarepickedupbyEEGvaryonamuchmorerapidtimescalethan thesluggishBOLDﬂuctuationsthat(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-speciﬁcHRF,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𝜽(𝑘)_into_a_full_HRF, representedasavector𝐡𝑘; isanoperatorthattransforms𝐡𝑘intoa Toeplitzmatrix𝐇𝑘bypopulatingthemainandlowerdiagonalswith theHRFsamples(seealsoA.2);𝑏 _𝑘,𝑖_𝑣istheweightforthe𝑘 -thHRFbasis functioninthe𝑖𝑣-thROI;𝐇𝑖𝑣 istheToeplitzmatrixholdingthetotal HRFinthe𝑖𝑣-thROI.ThisoperationisclariﬁedinFig.1b.

3_In_this_paper,_we_use_the_term_{‘neurovascular}_coupling’_to_describe_the cou-plingcharacteristicbetweenEEGandfMRItemporaldynamics,andmakethe silentassumptionthatthischaracteristicisaproxy/surrogatefor ‘neurovascu-larcoupling’asitisunderstoodinneuroscience:themodelthatdescribesBOLD changesinresponsetoelectricalneural‘events’,whichtaketheformoflocal ﬁeldpotentialsatthesynapses.

(5)

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𝜽(𝑘)_,_𝑘₌₁_…_𝐾 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𝑡,_where_we_omit_the_superscript(𝑘)_from_the pa-rameters𝜽 tonotoverloadthenotation.

2.5. Coupled matrix-tensor factorization of EEG and fMRI

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)

(6)

=[𝐇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 _HRF_basis _coeﬃcient _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-nonspeciﬁc relative weight for each basis func-tion𝑘(capturedin𝐛𝑘),andsource-speciﬁcamplitudes(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 coeﬃcients,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𝐽,composedoftwodataﬁttingtermsand 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 oftheresidualspromoteagoodﬁtofthelow-rankapproximationsto thedata,while the𝓁1-regularizationtermspenalizeexcessivesource amplitudesandpromoteaparsimonious4_model,_similar_to_the group-LASSOmethod(Acaretal.,2014;YuanandLin,2006).Atthesametime, thelatterpenaltyalsotendstopreventtheoccurrenceof degenerate terms(Bro,1997).Weminimize(14) usingtheStructuredDataFusion 4 _The_{sparsity-promoting}_properties_of_the_LASSO_penalty_are_most_useful_in thecontextofanunderdeterminedsystem,withmorecoeﬃcientsthan 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𝑅(seeAppendixB fordetails).

2.6. Statistical inference

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𝐿=250surrogatefMRĨ𝐘(𝑙)_datasets_using_the pro-cedureinBullmoreetal.(2001).Weresampleonlytheadjusteddata 𝐘−𝐍𝐏T,i.e.,afterremovingthecomponentswhich modelvariation specifictothefMRIdata.Weperforminferenceonapseudot-statistic, whichwecomputeforeveryROIandforeverysourceasfollows:

1. constructalocal‘designmatrix’withallestimatedtemporal signa-turesasin(3):𝐃𝑖𝑣 = [ 𝐳(1) 𝑖𝑣 …𝐳 (𝑅) 𝑖𝑣 ] , 2. ﬁndthenew‘betas’bysolving𝜷(𝑙)

𝑖𝑣 =𝐃 † 𝑖𝑣 ̃𝐲

(𝑙) 𝑖𝑣 ∀𝑙,

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

Throughthisprocedure,weobtain𝐿-pointempiricalnull distribu-tionsforeverysourceandeveryROI.Wesetthesigniﬁcancethresholdas tocontrolthefamilywiseerror(FWE)rateat 𝛼 = 0. 05, accordingtothe maximumstatisticprocedureoutlinedinNicholsandHayasaka(2003). Thatis,foreverysource𝑟,weformtheempiricaldistributionofthe max-imalt-statisticoverall𝐼 _𝑣ROIs,anddeterminesource-speciﬁcthresholds 𝑇(𝑟)

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

(𝛼)asthe 5%-percentile(totestfordeactivation).Finally,weobtainstatisticalmaps forallsources𝑟byapplyingthesethresholdstotheoriginalspatial sig-natures𝐯𝑟, whichcanbeconsideredasthebetasoftheunshuﬄeddata. Furthermore,wecreateamapoftheHRFvariabilityoverROIs.For everyROI,weassesshow‘unusual’thelocalHRFis,bymeasuringits calibrateddistanceinHRFspacetoallotherROIs’HRFs.Weusetwo metricstoquantifythis(seeAppendixCfordetailsonthecomputation). 1. Extremity iscomputedasoneminustheaverageoftheabsolute val-uesofthecorrelationsbetweenaHRFwaveformandallotherHRFs’ waveforms.

2. Entropy oftheHRFwaveformiscomputedasthenegativelogarithm oftheconditionalprobabilityoftheHRF.

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

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

2.7. Model performance

WeuseseveralmetricstoquantifythequalityoftheobtainedsCMTF solutions.

(7)

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-speciﬁc model selection

TableB.3 compilestheresultsofthemodelselectiondescribedin AppendixB.Foreachpatient,weselectthesetofsCMTFfactorsofrank ̂𝑅 , whichbestfulfillthecriteria.Inallcases,wefoundatleastonesucha solution,includinganIED-relatedcomponentwithinthatsolution.Note thatsometimesmodelswithdifferent𝑅mightscorewellondifferent (subsetsofthe)criteria,sotheselectionoftherankisinevitably am-biguous.Inthenextsection,weanalyzetheindividualsetofresultsfor

(8)

Fig.3.Goodnessofﬁtofeachpatient’sEEGtensor

andfMRImatrix𝐘,forvaryingchoicesoftherank𝑅

in thesCMTF. Naturally,theEEG approximation er-rordecreasesmonotonicallyforincreasingrank (intra-patient).ForthefMRIdata,theﬁtalreadyplateausfor verylow𝑅.Thisisduetothepresenceofadditional, uncoupledcomponents𝐧𝑞 ◦𝐩𝑞 inthefMRIfactorization, whichcanabsorbsomeofthevariancewhenthe num-berofcoupledcomponentsislow,butwhichlosetheir relevanceathigherranks.

eachpatient,basedontheselectedrank,andweanalyzethesensitivity oftheresultstothechoiceof𝑅 .

Weshowthegoodnessof ﬁtoftheestimatedfactorsfortheEEG tensor andthefMRImatrix𝐘inFig.3.Duetothenormalization stepswhichhavebeenappliedtothedata(cfr.Section2.2),thesCMTF operatesinaregimeofmoderatelyhighrelativeapproximationerrors. 3.2. Spatio-temporo-spectral proﬁles of interictal discharges

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).

Westressagainatthispointthatasubset oftheresultsisproneto 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,wealsoinvestigatethesigniﬁcantdeactivationsof thesources5_.

NotethatwedesignedtheHRFvariabilitymetricssothattheyareim- mune tothepolarityoftheHRFs.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. Patient 3

Weanalyzethesolutionwith ̂𝑅 =2sources,andshowtheresultsin Figs.4 and5.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.5 confirmsthis,and alsoshowsthattheIED-relatedsourceissignificantlyactiveindifferent ROIsoftheIOZ.

3.2.2. Patient 10

Weanalyzethesolutionwith ̂𝑅 =5sources,andshowtheresultsin Figs.6 and7.ThereisaclearIED-relatedsource,andalsoan artifac-tualsourceat±33Hz,whichisalsopresentinotherpatients.Duetoits relativelyconsistentoccurrence,wehypothesizethatthisartifactisdue totheMRacquisition.Forexample,itmaybearemnantofagradient

5_{Alternatively,} _it_is_possible_to_use_a _pseudo_F-statistic,_e.g. _the_squared pseudot-value,tobypassthesigncorrectionaltogether.Thedownsideofsuch anapproachisthatitisthenimpossibletodistinguishactivationand deactiva-tion,whichmaybemeaningful.

(9)

Fig.4. (a)Intheselectedsolutionforpatient3(̂𝑅=2),bothsourceshaveatemporalsignaturethatcorrelatedstronglytothereferenceIEDtimecourse.Theﬁrst sourcemodeledthemainonsetofIEDsandwaslow-frequencyandtopographicallyfocal,whilethesecondsourcewasspatiallyandspectrallydiﬀuseandcaptured thepropagationofIEDstoremoteareas(cfr.Fig.8).(b)FiveoutofthetwentymostdeviantHRFswerefoundinsidetheictalonsetzone(boldlines,𝑝<10−4_)._These HRFshadmainpeaksbefore0s,i.e.,theyledtoBOLDchangesbeforethecorrespondingEEGcorrelateoftheIEDwasseen.

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

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

3.2.3. Summary of all patient’s results

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

(10)

Fig.5.ThestatisticalnonparametricmapsoftheIED-relatedcomponent(toptworows)andHRFentropy/extremitymaps(bottomtworows)ofpatient3show concordancewiththeictalonsetzone(IOZ).EspeciallytheregionsofsigniﬁcantIEDactivationwereaccurate,butalsoﬁveoutofthetwentyregionswiththemost deviant(highestentropy)HRFswerefoundintheIOZ(cfr.Fig.4b).Thegroundtruthictalonsetzoneishighlightedindarkgraywithawhitecontour.ROIswith highvaluesforbothHRFvariabilitymetricsarecoloredinorange.

Table2

ThesCMTFleadstothreetypesofspatialinformation,whichcanbecross-validatedagainstthegroundtruthIOZ,asdeﬁnedinSection2.7andsummarizedforall patientsinTable1:(1)theEEGtopography𝐦IEDoftheIED-relatedcomponent;(2)thesigniﬁcantlyactivatedanddeactivatedROIsinthefMRIspatialsignature

𝐯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-ablebutstructuredﬂuctuationsathigherfrequenciesmaybecaptured inacomponent.

3.2.4. Sensitivity to model selection

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.

(11)

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

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

A novel EEG–fMRI data fusion framework

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-speciﬁcpreprocessingandtensorization steps,whichleadtoathird-orderEEGspectrogramtensorvaryingover electrodes,timepoints,andfrequencies,andanfMRImatrixwithBOLD timecoursesforapredeﬁnedsetofregionsofinterestorparcels;(3) coupledmatrix-tensorfactorizationoftheEEGtensorandfMRImatrix alongthesharedtemporalmode,whileaccountingforvariationsinthe localneurovascularcoupling;(4)automatedselectionofarobust,and relevantIED-relatedcomponent,andnonparametrictestingtoinferits spatialdistributioninthebrain.

Wehavestressedtheimportanceofandaccountedforthevariability of thehemodynamicresponsefunction (HRF)over diﬀerentpatients

(12)

Fig.7.ThestatisticalnonparametricmapsoftheIED-relatedcomponentandHRFentropy/extremitymapsofpatient10showconcordancewiththeictalonsetzone (IOZ).IEDoccurrenceswereassociatedwithsigniﬁcantlyactive(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 havecomputeddiﬀerentsummarymetrics(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, wealsofoundregionswithintheIOZthatsigniﬁcantlydeactivatedin 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,usingdiﬀerenttypesof EEG-derived

(13)

regres-S. Van Eyndhoven, P. Dupont, regres-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 oftendistinctfromtheROIsidentiﬁedassigniﬁcantlyactivatedtothe 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.

HRFs vary strongly over subjects and brain regions

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,itstrictlymeansthatwe observed BOLD changesinthefMRIdatathatoccurbeforethecorresponding observed neuralcorrelateontheEEG.DespitethefactthatmanyoftheHRFs dif-feredsubstantiallyfromthecanonicalHRF,whichiscausalandpeaks approximately 6s after 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.

Prior EEG signal enhancement aids analysis

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).

The interpretation of components

Overall,thesCMTFpipelinesucceededinextractingmeaningful IED-relatedcomponents,alongsidecomponentsthatmodeled resting-state neuralﬂuctuationsandphysiologicalandtechnicalartifacts.Thefact thatthesCMTFcanestimatesignaturesandstatisticalmapsformultiple componentsisapowerfuladvantageoverclassicalEEG-correlated anal-ysis.Aswedemonstratedintheexperiments,artifactualinﬂuencesmay