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
Huffel
a,
Borbála
Hunyadi
ga 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 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 wereacquired on one of two 3T MR scan-ners(AchievaTXwitha32-channelheadcoilandInteraAchievawith aneight-channelheadcoil,PhilipsMedicalSystems,Best,The Nether-lands)withanechotime(TE)of33ms,arepetitiontime(TR)of ei-ther 2.2or2.5s,andavoxel sizeof2.6×3×2.6mm3 . EEGdata 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).TheseprocessingstepswerecarriedoutusingSPM8 (Func-tionalImaging Laboratory,Wellcome Centerfor Human Neuroimag-ing,UniversityCollegeLondon,UK)(Fristonetal.,1994).Wereferthe readertoTousseynetal.(2014a)foradetaileddescriptionofthese pre-processingsteps.
WeregressoutcovariatesofnointerestfromthefMRIdata.These include:thesixmotion-correctionparameters(plustheirsquaresand 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 data of thedifferent runs aretemporally concatenated. Further cus-tomizedpreprocessingstepsaretreatedinSections2.3 and2.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)or EEG 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 arethenfiltered as𝐖̂T ̃𝐱.Duetotheextensionwithlaggedcopiesof thesignals,channel-specificfiniteimpulseresponsefilters arefound. Hence,𝐖̂T ̃𝐱isasetofspatiotemporallyfilteredoutputsignals,inwhich IED-like waveformsarepreservedwhile otherwaveforms,which are notspecifictoepilepsy,aresupressed1 .
WetraintheMWFforeach patientindividually,afterembedding theEEGsignalsusing𝜏 =4positiveandnegativelags2 ,andcompute thefinalfilterusingthegeneralizedeigenvaluedecomposition,which ensuresthepositivedefinitenesspropertyofthesubtractedcovariance matrixin(2) (Somersetal.,2018).
1Toretrievefilteredversionsoftheoriginalsetofchannelsonly(andnotof theirtime-delayembeddedcopies,whichwouldberedundant),onlythemiddle
𝐼𝑚 columnsof𝐖̂ areused(cfr.(1)).
2WeobservedinVanEyndhovenetal.(2019a)thatformostpatients,the outputSNRsaturatesaroundthisvalue,correspondingtoanintervalof-16ms to+16ms.
Table1
Clinicalpatientdata.Abbreviations:F=female,M=male,L=left,R=right,CNS=centralnervoussystem,DNET=dysembryoplasticneuroepihelial tumor,FCD=focalcorticaldysplasia,HS=hippocampalsclerosis,cx=cortex,IED=interictalepilepticdischarge,TR=repetitiontime,IEDloc.= localizationoftheIEDonEEG.
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
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.
2.4.1. Spatio-temporal-spectraltensorrepresentationofEEG
Weadopta tensorizationstrategybasedon time-frequency trans-formationoftheEEGdatatothird-orderspectrograms(timepoints× frequencies×channels).Afterthepre-enhancementstepdescribedin Section2.3, we createa spectrogram using theThomsonmultitaper 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)
3Inthispaper,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-orderspectrogramtensorisfactorizedaccordingto(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𝐯𝑟 .
Here,𝐬𝑟isafactorvectorholdingthetimecourseofthe𝑟-thEEGsource; isanoperatorthattransformsasetofparameters𝜽( 𝑘) intoafullHRF, representedasavector𝐡𝑘; isanoperatorthattransforms𝐡𝑘intoa Toeplitzmatrix𝐇𝑘 bypopulatingthemainandlowerdiagonalswith theHRFsamples(seealsoA.2);𝑏𝑘,𝑖𝑣 istheweightforthe𝑘-thHRFbasis functioninthe𝑖𝑣-thROI;𝐇𝑖𝑣 istheToeplitzmatrixholdingthetotal HRFinthe𝑖𝑣-thROI.ThisoperationisclarifiedinFig.1b.
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⋅ 𝜃𝜃1
2 𝑡𝜃1−1𝑒−𝜃2𝑡−
𝜃5 Γ(𝜃3 )−1 ⋅ 𝜃𝜃3
4 𝑡𝜃3−1 𝑒− 𝜃4𝑡,whereweomitthesuperscript( 𝑘) fromthe pa-rameters𝜽 tonotoverloadthenotation.
2.5. Coupledmatrix-tensorfactorizationofEEGandfMRI
Aftertensorization,wejointlydecomposetheEEGtensorandthe 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(alsocalled‘component’) intermsof an outerproduct(◦)ofatemporal,spectral,andspatialsignature, respec-tively.Unlikematrixdecompositions, thedecomposition ofa 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) wherê 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 𝐹 +𝛾𝑥||𝝀𝑥||1 + 𝛾𝑦||𝝀𝑦||1 (14) s.t.𝐇𝑘=(𝐡𝑘)=((𝜽( 𝑘) )) 𝝀𝑥=[𝜆𝑥,1 … 𝜆𝑥,𝑅]T 𝜆𝑥,𝑟=||𝐬𝑟||2 ⋅ ||𝐠𝑟||2 ⋅ ||𝐦𝑟||2 𝝀𝑦=[𝜆𝑦,1 … 𝜆𝑦,𝑅]T 𝜆𝑦,𝑟=∑𝐾𝑘=1 ||𝐛𝑘∗𝐯𝑟||2 , (15)
wherethesquaredFrobeniusnorm||||𝐹2 ofatensoristhesumof itssquaredelements;||𝐚||2 and||𝐚||1 denotetheEuclideanor𝓁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, thelatterpenalty alsotendstopreventtheoccurrenceofdegenerate terms(Bro,1997).Weminimize(14)usingtheStructuredDataFusion frameworkinTensorlab(Sorberetal.,2015;Vervlietetal.,2016),using aquasi-Newtonmethodbasedonalimited-memoryBFGSalgorithm,for 50independentinitializations(seeAppendixAfordetailsregardingthe optimizationprocedureandparameters).Afterconvergence,eachsetof estimatedfactorsneedstobecalibratedtoremovecertainambiguities, andmodelselectionmustbeperformedtopickthebestsolution,with anappropriate𝑅(seeAppendixB fordetails).
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.
4Thesparsity-promotingpropertiesoftheLASSOpenaltyaremostusefulin thecontextofanunderdeterminedsystem,withmorecoefficientsthan observa-tions,e.g.indictionarylearning.Here,theproblemisheavilyoverdetermined, andwedonotexpectthattheamplitudes𝝀𝑥 an𝝀𝑦 goexactlytozero.However, the𝓁1-penaltyisstillausefulheuristictoavoiddegeneratecomponentsinthe EEG’sCPdecomposition.
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.
2. EntropyoftheHRFwaveformiscomputedasthenegativelogarithm oftheconditionalprobabilityoftheHRF.
Bothforthepseudot-mapsasfortheHRFextremityandentropy maps,wefurthermorelimittheinspectiontothe20ROIswiththe high-estvalues,ifapplicable.
Anend-to-endoverviewofourpipeline,fromdatapreprocessingup untilstatisticalinference,isdepictedinFig.2.
2.7. Modelperformance
WeuseseveralmetricstoquantifythequalityoftheobtainedsCMTF solutions.
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 voxelsorregionsoutsideof thedelineatedIOZasfalsepositives. Ac-knowledgingepilepsyasanetworkdisorder,suchactiveregionsmight reflectseizureorIEDpropagation,despitenotbeinginvolvedintheir generation.
Furthermore,wehypothesizethatthespatialvariationoftheHRF over thebrainmightreveal additionallocalizinginformation regard-ingtheIOZ,i.e.,basedonconsiderationsexplained inSection1,we assumethattheHRFinorneartheIOZmightbedistortedcomparedto nonepilepticbrainregions.Wetestthishypothesisbyassessingwhether thoseregionscorrespondtohighvaluesintheHRFentropyandHRF ex-tremitymaps(cfr.Section2.6).
Additionally,weinspectthespectral,spatialandtemporalEEG sig-naturesoftheextractedsources,andwemeasurewhetherthespatial fMRIsignaturesbearanysimilaritytoknownnetworksofresting-state humanbrainactivity(Shireretal.,2012).
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.
3. Experiments
3.1. Patient-specificmodelselection
TableB.3 compilestheresultsofthemodelselectiondescribedin AppendixB.Foreachpatient,weselectthesetofsCMTFfactorsofrank
̂𝑅,whichbestfulfillthecriteria.Inallcases,wefoundatleastonesucha solution,includinganIED-relatedcomponentwithinthatsolution.Note thatsometimesmodelswithdifferent𝑅mightscorewellondifferent (subsetsofthe)criteria,sotheselectionoftherankisinevitably am-biguous.Inthenextsection,weanalyzetheindividualsetofresultsfor 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𝐬ref forselectingtheIED-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,thisis morechallengingforfMRI.Thatis,frequently,the polarityoftheHRFwaveformisambiguous,andmakingthe‘wrong’ choicein avoxel𝑖𝑣 (i.e.,theHRFhastheoppositeeffectofthetrue physicalcerebralbloodflowchange)immediatelyleadstothewrong signof thespatialcoefficientsin𝐯𝑟intherespectivevoxel,andtheir pseudot-values,forallsources𝑟.Totracktheoccurrenceofthis fore-seenfailuremode,wealsoinvestigatethesignificantdeactivationsof thesources.5
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.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,
5Alternatively, 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.
i.e.,withapositivepeakbeforezeroseconds.Fig.5 confirmsthis,and alsoshowsthattheIED-relatedsourceissignificantlyactiveindifferent ROIsoftheIOZ.
3.2.2. Patient10
Weanalyzethesolutionwith ̂𝑅=5sources,andshowtheresultsin Figs.6 and7.ThereisaclearIED-relatedsource,andalsoan artifac-tualsourceat±33Hz,whichisalsopresentinotherpatients.Duetoits relativelyconsistentoccurrence,wehypothesizethatthisartifactisdue
totheMRacquisition.Forexample,itmaybearemnantofagradient artifactwhichisnotadequatelyremovedfromthedataofsome chan-nels,cfr.theobservationmadeinMarečeketal.(2016).Surprisingly, thissourceissignificantlyactivein anextendedareain theoccipital lobe,overlappingwiththevisualnetwork.BothHRFmetricsreached extremevalues atsome(distinct)ROIs withintheIOZ.Thepseudo-t mapoftheIED-relatedsourceshowssignificantlyactiveROIsthatare concordantwiththeIOZ,anddeactivationofalargepartofthedefault modenetwork.Furthermore,theIED-relatedsource’sEEGtopography
Fig.5.ThestatisticalnonparametricmapsoftheIED-relatedcomponent(toptworows)andHRFentropy/extremitymaps(bottomtworows)ofpatient3show concordancewiththeictalonsetzone(IOZ).EspeciallytheregionsofsignificantIEDactivationwereaccurate,butalsofiveoutofthetwentyregionswiththemost deviant(highestentropy)HRFswerefoundintheIOZ(cfr.Fig.4b).Thegroundtruthictalonsetzoneishighlightedindarkgraywithawhitecontour.ROIswith highvaluesforbothHRFvariabilitymetricsarecoloredinorange.
isveryconsistentwiththeclinicaldiagnosis.Thefourthsourceisactive inthedefaultmodenetwork,predominantlyinthe𝛼 band(cfr.Fig.9). Thefifthsourcehadanunclearspectrum,butitstemporalsignature correspondstotheoccurrenceofhigh-amplitudeIEDs.Itspseudot-map showswidespreadactivationsoverthebrain,whichdidnotincludethe IOZ.WeexpectthatthiscomponentcapturesthepropagationofIEDs, afteronsetneartheIOZ,similarlytopatient3.
3.2.3. Summaryofallpatient’sresults
Weprovideanoverviewoftheresultsw.r.t.IOZdetectioninTable2. Allresultstakentogether,thesCMTFresultsallowacorrectdetection oftheIOZbasedonthesignificantIEDactivation(10/12cases)and 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 whichis nextin line,i.e.,which wouldbeasecondbest (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.
4. Discussion
AnovelEEG–fMRIdatafusionframework
We have proposed an integrated andstructured coupled matrix-tensorfactorization(sCMTF)framework,whichcan beusedtomake 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 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-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). 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.Inroughlyhalf ofallcases, 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-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). tionsofthemodel.Hence,theoutputofthisanalysisisonlytobeused
inconjunctionwithothermodalities(e.g.SISCOM)duringpresurgical evaluation,inordertoassesscross-modalityagreement,asis already commonforEEG–fMRI.
Wepreviouslystudiedthesamepatientcohortusingclassical EEG-correlatedfMRIanalysis,usingdifferenttypesofEEG-derived regres-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 11 hadonlyfewIEDs(15and6,respectively),which makesmining of theirEEG–fMRI dataforBOLDactivationdifficultapriori( Salek-Haddadietal.,2006;2003;Zijlmansetal.,2007).Hence,alsoforour moreflexiblemethod,theprobableyieldscaleswiththenumberofIEDs duringrecording.Weinspectedthe20HRFsandROIswiththe high-estextremityandentropy.Hence,itisinevitablethatsomeormostof theseROIsarenotwithintheIOZ,ortheIOZmightevennothavea deviant HRF.StandaloneHRFmetricswouldhencehaveahighfalse discoveryrate,eventhoughforseveralpatients,thehighproportionof IOZ-coveringROIsamongthe20selectedROIswasveryunlikelydueto chance(asmeasuredwithp-values).Still,itisnotalwaysthecasethat theHRFinsidetheIOZismeasurablydifferentthanintherestofthe brain,inwhichcasethesemetricswouldnotbeverysensitive.
However,theROIsthatwerehighlightedbytheHRFmetricswere oftendistinctfromtheROIsidentifiedassignificantlyactivatedtothe
S. Van Eyndhoven, P. Dupont, S. Tousseyn et al. NeuroImage 228 (2021) 117652 Table 2 The sCMTF leads to three types of spatial information, which can be cross-validated against the ground truth IOZ, as defined in Section 2.7 and summarized for all patients in Table 1 : (1) the EEG topography 𝐦IED of the IED-related component; (2) the significantly activated and deactivated ROIs in the fMRI spatial signature 𝐯IED ; (3) the ROIs with strongly deviating HRF waveforms, as measured via entropy and extremity. Since the EEG topography has a very low spatial resolution, and depends on the attenuation properties of the tissue as well as the orientation of the neural sources in the cortex, we only expect partial similarity to the IOZ’s spatial focus; hence, we use the term ‘consistent’ rather than ‘concordant’. The patients who had a good outcome after surgery (patients 5, 7 and 11) had a higher concordance between the three types of spatial clues than patients with a poorer outcome (patients 1, 2 and 6). patient select e d solution EEG t o pogr aph y consis te nt with IOZ? signatur e 𝐯 IED concor d ant with IOZ? HRF va ri a b il it y me trics complement ar y to 𝐯IED ? 20 highes t-entr op y RO Is ID ̂ 𝑅 acti v ation deacti v a tion entr op y ex tr e m it y # in IOZ (p-v alue) 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 )
IEDs.Hence,theSnPMsoftheIEDsandtheentropymetricsprovidevery complementaryinformation,andwhenanalyzedjointly,theymayinfer thelocationoftheIOZwithmore certainty,byselectingbrainareas 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-rent methodcouldthenallowtoextractRSNswhicharereflectedin bothEEGandfMRIdata,giventhatsufficientlymanycomponents𝑅
areextracted.
HRFsvarystronglyoversubjectsandbrainregions
Thereweresubstantialdifferencesin(estimated)neurovascular cou-plingoverpatientsandbrainregions,asexpected.Sinceweused ‘regu-larized’basisfunctions,whichareparametrizedassmoothgamma den-sityfunctions,theresultingHRFsgenerallyhadaplausibleshape. How-ever, insomecaseswe foundnonsensicalshapes,in which,e.g., the waveformahadthesamepolarityoverthewholetimecourse, poten-tially withabimodalshape (cfr.patient4).Thisservesasa humble 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.,datafrom thesamepatient,recordedoutside theMR scanner). WhilethistaskstillfrequentlyreliesontheskillofhumanEEGreaders andneurologists,advancedautomatedsolutionsforinterictalspike de-tectionareavailable(Scheueretal.,2017;Wilsonetal.,1999).Within eachsolutionofaspecificrank,wepickedtheIEDcomponentasthe onewiththehighestcorrelationwithareferencetimecoursedirectly derivedfromtheenhancedEEG.Someofthepresentedresultsmake clearthatthisreferencetimecourseis notcompletelyfreefrom 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,alongsidecomponentsthatmodeledresting-state neuralfluctuationsandphysiologicalandtechnicalartifacts.Thefact thatthesCMTFcanestimatesignaturesandstatisticalmapsformultiple componentsisapowerfuladvantageoverclassicalEEG-correlated anal-ysis.Aswedemonstratedintheexperiments,artifactualinfluencesmay beisolatedinseparatecomponents,whichcouldreduce theirimpact onIEDmappinginthebrain.Additionally,weencounteredcaseswhere twocomponentswerecorrelatedtotheIEDoccurrences:thecomponent withthehighesttemporalcorrelationtoareferenceIEDtimecoursethen correctlyrevealedthelocalizationoftheIOZ,whiletheothercomponent presumablymodeledthepropagationofIEDstoremotebrainregions. ThisobservationisanalogoustothefindinginHunyadietal.(2016), whereadifferenttypeofCMTFwasappliedtoaverageEEGwaveforms ofIEDsandstatisticalBOLDmaps,whichrevealedadissociation be-tweentheearlyIEDspikeandthesubsequentwave,whichwererelated totheonsetandspreadoftheIEDs,respectively.Sincewetransformed thedatawithatime-frequencytransformthatusedwindowsoflength TR,ouralgorithmisunabletounraveldifferentphaseswithinoneIED, sincetheyoccurinashortertimeframe.However,weidentifiedthese differentIED-relatedsourcesbytheirsignificantlycorrelatedtemporal signatures,andtheirdistinctspatialandspectralprofiles.Whilewedid notimposenonnegativityconstraints,manyestimatedEEGspectraland spatialsignatureswereapproximatelynonnegative.Thisneednotbethe case,however,sincetheEEGdataarenormalizedinawaysuchthatthe resultingsignatureswouldrevealrelativeincrease/decrease,ratherthan absolutetime-varyingspectralpower(Marečeketal.,2016).For exam-ple,ifacertaincomponentisassociatedwithapowerincreaseinone spectralband,andasimultaneouspowerdecreaseinanotherband,this wouldbereflectedinaspectrumwithbothapositiveandanegative peak.
Practicalconsiderations
Theend-to-endsCMTFpipelinecanprovidearichersetofresults comparedtoclassicalEEG-correlatedfMRIanalysis.Inthisrespect,itis amorepowerfuldataexplorationtool.Thetradeoff tobemadeisthat
significantcomputationtimegoesintothesCMTFandsubsequent in-ference—ifonewantstoapplyitasrigorouslyaswehavedoneinthe currentexperiments.Weseemtobedoingalotofunnecessarywork, bycomputingthesCMTFfactorsforseveralnumbersofsources,andby repeatingtheoptimizationseveraltimesforafixednumberofsources. Unfortunately,bothwaysofrepetitionseemrequiredtoobtainrobust results,aswehavearguedinAppendixAandAppendixB.However,the EEG-onlyCPdecomposition,whichliesattheheartofourinitialization strategy,seemedveryrobust:wefoundhighlysimilarEEGsignatures foralmostallrandominitializations.Probably,thisisthankstotheuse ofthepowerfulGauss–Newton-typeoptimization.Hence,fewer repeti-tionsofthesCMTFmaybealreadysufficienttoarriveatthesamerobust results.DespitetheveryreproducibleEEGsignaturesintheinitialCP decompositions,westillperformed50repetitionsofthesCTMF,each timeslightlyvaryingtheinitialHRFparameters.Assuch,webelieve ourfindingsarereasonablyrobusttopoorinitializationoftheHRFs. Per-formingthesCMTFformanychoicesof𝑅maystillberequired,asthe qualityoftheresultdependsontheextractionofanappropriate num-berofsources.Inourstudy,noprohibitivecomputationswereneeded, sincetheheuristicselectionprocedurepreferredlowranks,whichwas stillsufficienttomodeltheIED-relateddynamics.Toalsocapturemore resting-stateactivity,thecandidateranksandthemodelselection pro-cesscouldbechosendifferently.Furthermore,wehavedemonstratedin ourexperimentsthatthesummarymetrics(sensitivityforlocalizingthe IOZbasedondifferentstatisticalscores)arefairlyrobusttothechoice of𝑅,althoughtheestimatedsignaturesthemselvesdiffer.
Formanypatients,theavailabledataweresplitacrossmultipleruns (i.e.,withafewminutesbreakinbetween),andweoptedto tempo-rallyconcatenatedataoverruns,asexplainedinSection2.2.Whilethis violatesthecouplingmodelbasedonHRFsfortimesamplesnearthe boundaries,weconsidertheeffectminimal,giventhatthenumberof those‘affected’timesamplesrepresentsaverytinyfractionofthewhole timeseries.However,amorerigorousapproachwouldbeto ‘inverse-impute’ thesesamplesandconsiderthemasmissingvalues:assuch, theyareignoredduringthesCMTFoptimizationandwillnotaffectthe results.
Strategiestoalleviatethecomputationaldemand
Duetotherepeateddecompositionandthenatureofthe nonpara-metricinference,thecomputationsarehighlyparallellizable.Fora typ-icaldatasetwithavailableIEDannotations,andwiththeparameterswe haveused forthisstudy,theend-to-endcomputationforonepatient tookatmostfivehoursonamachinewithtwelvecores.Toalleviatethe computationalburden,wehaveparcellatedthefMRIdatainto246 re-gions,basedontheBrainnetomeatlas(Fanetal.,2016).Thisisclearly suboptimal,astheatlasisnotpatient-specific,andismostlydesignedto studyhealthybrains.Thereisaseriousriskforpartialvolumeeffects, inwhichtheIOZisscatteredoverseveralROIs.Assuch,theIED-related BOLDchangesinthepartoftheIOZthatfallswithinacertainROImay getswampedbytheremainingBOLDfluctuationswithintheROI de-lineation.Hence,wehopetobeableovercomethisproblem,eitherby algorithmicimprovements,includingaspeedupoftheoptimization,or bytheuseofapatient-specificparcellationorPCA-likecompressionof thefMRIdata.Asofyet,itishardtosaywhetherthefixedatlashadan adverseeffectontheresults,anditisnotsostraightforwardtocompare thestatisticalmapsfromthisstudytomapswhicharevoxel-based.We arecurrentlypursuingexperimentsinwhichweemployahierarchical parcellation:inafirststep,theBOLDtimeseriesaregrouped(butnot yetaveraged)accordingtotheBrainnetomeatlas;subsequently,weuse spectralclusteringtofurtherrefineeachBrainnetomeparcelbasedon thecorrelationmatrixofitsBOLDtimeseries.Assuch,thishybrid ap-proachcombinesafixed,coarse-grainedatlaswithafurtherdata-driven subdivision,whichcanmitigatepartialvolumeeffects,whilestill pro-vidingasignificantdatacompression.Forpatientswithlesions,a cus-tomized parcellationcan beused,inwhichthelesionitselfcoincides withoneparcelorwiththeunionofseveralparcels.Alternatively,it
