structuredcmtf
Pipeline for structured coupled matrixtensor factorization of EEGfMRI data
view repo
EEGcorrelated fMRI analysis is widely used to detect regional blood oxygen level dependent fluctuations that are significantly synchronized to interictal epileptic discharges, which can provide evidence for localizing the ictal onset zone. However, such an asymmetrical, massunivariate approach cannot capture the inherent, higher order structure in the EEG data, nor multivariate relations in the fMRI data, and it is nontrivial to accurately handle varying neurovascular coupling over patients and brain regions. We aim to overcome these drawbacks in a datadriven manner by means of a novel structured matrixtensor factorization: the singlesubject EEG data (represented as a thirdorder spectrogram tensor) and fMRI data (represented as a spatiotemporal BOLD signal matrix) are jointly decomposed into a superposition of several sources, characterized by spacetimefrequency profiles. In the shared temporal mode, Toeplitzstructured factors account for a spatially specific, neurovascular `bridge' between the EEG and fMRI temporal fluctuations, capturing the hemodynamic response's variability over brain regions. We show that the extracted source signatures provide a sensitive localization of the ictal onset zone, and, moreover, that complementary localizing information can be derived from the spatial variation of the hemodynamic response. Hence, this multivariate, multimodal factorization provides two useful sets of EEGfMRI biomarkers, which can inform the presurgical evaluation of epilepsy. We make all code required to perform the computations available.
READ FULL TEXT VIEW PDFPipeline for structured coupled matrixtensor factorization of EEGfMRI data
Refractory epilepsy is a neurological disorder suffered by 30% of approximately 50 million epilepsy patients worldwide (who2019epilepsy2), in which seizures cannot adequately be controlled by antiepileptic medication. In the preparation of treatment via resective surgery, interictal epileptic discharges (IEDs) can be localized in the brain with simultaneous EEGfMRI, which provides a good surrogate for mapping the seizure onset zone (lemieux2001event; thornton2010eeg; vanhoudt2013eeg; grouiller2011or; zijlmans2007eeg; an2013electroencephalography; khoo2017hemodynamic). This mapping is often conducted via EEGcorrelated fMRI analysis, wherein a reference temporal representation of the IEDs is used to interrogate all brain regions’ blood oxygen level dependent (BOLD) signals for significant correlations; voxels for which a statistical threshold is exceeded can then be considered part of the epileptic brain network, along which epileptic seizures are generated and propagated (gotman2008epileptic; lemieux2001event; zijlmans2007eeg; thornton2010eeg; salek2003studying).
The workhorse for conducting EEGcorrelated fMRI analyis has been (salek2006hemodynamic)—and will likely continue to be (poline2012general)—the general linear model (GLM) framework (friston1994statistical)
. Over the past years, it has become clear that using the GLM comes with several hurdles, related to the many modeling assumptions, that may reduce its sensitivity or specificity (increasing Type I errors) when violated
(poline2012general; monti2011statistical; lindquist2009modeling). Remedies for several of these issues are not yet widely applied, or are not yet available.First of all, the adoption of a relevant representation of IED occurrences to construct a regressor for the design matrix has proven vital to the sensitivity. This aspect has been investigated in (rosa2010estimating; murta2015electrophysiological; abreu2018eeg; vaneyndhoven2019semi). In previous work (vaneyndhoven2019semi)
, we addressed this issue by preenhancing the EEG signals using a spatiotemporal filter that is tuned to maximize the signaltonoise ratio (SNR) of IEDs with respect to the background EEG. We have shown that taking the timevarying power of the filtered EEG leads to a robust regressor, which is more performant than many other types of regressors, including those based on stick functions
(lemieux2001event; salek2006hemodynamic), ICA (formaggio2011integrating; abreu2016objective) or EEG synchronization (abreu2018eeg).Model mismatch may occur due to the unknown neurovascular coupling from electrophysiological phenomena measured on the EEG to hemodynamic variations captured by the BOLD signals. In many papers on EEGcorrelated fMRI, a canonical hemodynamic response function (HRF) based on two gamma density functions is used to translate IEDrelated temporal dynamics to BOLD fluctuations (friston1998event). However, there is insurmountable evidence that the HRF is not fixed, but varies substantially over subjects (aguirre1998variability), over brain regions (handwerker2004variation), with age (jacobs2008variability), or even with stress level (elbau2018brain). For the diseased brain, this issue may be even greater: i.e., additional variation, e.g. in brain areas involved in the epileptic network, has been observed compared to healthy controls (vanhoudt2013eeg; benar2002bold; jacobs2009hemodynamic; lemieux2008noncanonical; grouiller2010characterization)
. Plenty of previous research has shown that failing to account for this variability may lead to substantial bias and increased variance of the estimated activation, which in turn inflates Type I and/or Type II error rates
(lindquist2009modeling; lindquist2007validity; calhoun2004fmri; monti2011statistical).Several methods have been devised to deal with this variability. A widely used approach is to model the HRF as a linear combination of several basis functions. Some popular choices for these bases, which are also supported by open source toolboxes like SPM are the ‘informed basis set’
(friston1998event), consisting of the HRF plus its derivative w.r.t. time and its derivative w.r.t. the dispersion parameter (leading to a Taylorlike extension which can capture slight changes in peak onset and width), and the finite impulse reponse (FIR) basis set, in which every basis function fits exactly one sample of the HRF in every voxel (glover1999deconvolution; aguirre1998variability). Other researchers have aimed to find a basis set by computing a lowdimensional subspace of a large set of ‘reasonable’ HRFs (woolrich2004constrained) or by fitting nonlinear functions to given fMRI data (lindquist2007validity; vaneyndhoven2017flexible). Alternatively, multiple copies of a standard HRF, which differ only in their peak latencies, can be used (bagshaw2004eeg). Finally, approaches exist that aim to be immune to differences in neurovascular coupling, such as those based on mutual information (MI), which does not rely on any predefined model or even linearity of the HRF (ostwald2011information; caballero2013mapping). Perhaps surprisingly, the authors of (caballero2013mapping) found that the results based on MI were often very similar to those based on the informed basis set, leading to the conclusion that the assumption of a linear timeinvariant system, as described by the convolution with an appropriate HRF, is sufficiently accurate. Instead, it may be useful to not make abstraction of the variable neurovascular coupling, but rather consider it as an additional biomarker to localize epileptogenic zones (vanhoudt2013eeg). Indeed, in several studies HRFs that deviate from the canonical model were found in regions of the epileptic network (benar2002bold; lemieux2008noncanonical; hawco2007bold; pittau2011changes; jacobs2009hemodynamic; moeller2008changes; vanhoudt2013eeg). Several hypotheses have been postulated to explain this varability, including altered autoregulation due to higher metabolic demand following (inter)ictal events (schwartz2007neurovascular), vascular reorganization near the epileptogenic region (rigau2007angiogenesis), or the existence of prespike changes in neuroelectrical activity which are not visible on EEG and which culminate in the IED (jacobs2009hemodynamic). It is thus an opportunity to map not only regions with statistically significant BOLD changes in response to IEDs, but also the spatial modulation of the HRF itself, in order to discover regions where an affected HRF shape may provide additional evidence towards the epileptic onset.The previous considerations indicate that it is difficult to meet all assumptions in the general linear model, which may compromise inference power (lindquist2009modeling; handwerker2004variation; monti2011statistical). Datadriven alternatives may relieve this burden, since they adapt to the complexity of the data more easily compared to modelbased approaches, and are especially suited for exploratory analyses (mantini2007electrophysiological; marecek2016can). Blind Source Separation (BSS) techniques consider EEG and/or fMRI data to be a superposition of several ‘sources’ of physiological activity and nonphysiological influences. Based on the observed data alone, BSS techniques are used to estimate both the sources and the mixing system, by means of a factorization of the data into two (or more) factor matrices, holding sources or mixing profiles along the columns. They naturally allow a symmetrical treatment of EEG and fMRI data, enabling true fusion of both modalities (valdes2009model; lahat2015multimodal; calhoun2009review), which is in contrast to EEGcorrelated fMRI, where EEGderived IEDs inform the fMRI analysis. Furthermore, BSS techniques naturally accommodate higherorder representations of the data in the form of tensors or multiway arrays, which can capture the rich structure in the data. Indeed, measurements of brain activity inherently vary along several modes (subjects, EEG channels, frequency, time, …), which cannot be represented using matrixbased techniques like ICA without loss of structure or information (sidiropoulos2017tensor; lahat2015multimodal; kolda2009tensor; acar2007multiway). Tensorbased BSS techniques have been used to mine unimodal EEG data by decomposing thirdorder spectrograms (channels time points wavelet scales) into several ‘atoms’ (also coined ‘components’ or ‘sources’), each with a distinct spatial, temporal and spectral profile/signature (miwakeichi2004decomposing; morup2006parallel; marecek2016can), with successful application in seizure EEG analysis (acar2007multiway; devos2007canonical). While a tensor extension of ICA for group fMRI data (in the form of subjects time points voxels) exists (beckmann2005tensorial), matrix representations of fMRI remain dominant for singlesubject analyses. Coupled BSS techniques can estimate components which are shared between both modalities, providing a characterization in both domains (hunyadi2017tensor). For example, in (acar2017acmtf; acar2019unraveling; hunyadi2016fusion; chatzichristos2018fusion), multisubject EEG and fMRI data have been analyzed using coupled matrixtensor factorization (CMTF), wherein the ‘subjects’ factor is shared between the EEG trilinear tensor decomposition and the fMRI matrix decomposition. In (hunyadi2016fusion), the resulting factor signatures revealed onset and propagation zones of an interictal epileptic network that was common over patients, as well as the modulation of the defaultmode network (DMN) activity. Also singlesubject data can be decomposed into distinct components, using a shared temporal factor for EEG and fMRI. This requires the use of a model of the neurovascular coupling, to ensure temporal alignment of EEG and BOLD dynamics. In (martinez2004concurrent), a fixed canonical HRF was used, followed by multiway partial least squares to extract components with spatial, temporal, and spectral signatures. In previous work, we proposed an extension to this technique, where a subjectspecific HRF is coestimated from the available data, along with the components (vaneyndhoven2017flexible).
In this paper, we extend this latter technique in order to account not only for subjectwise variation of the HRF, but also capture variations over brain regions. This results in a highly structured CMTF (sCMTF) of the interictal multimodal data, in which HRF basis functions and spatial weighting coefficients are estimated along with spatial, spectral and temporal signatures of components. By preprocessing the EEG using the datadriven filters from (vaneyndhoven2019semi), we aim to maximize the sensitivity in mapping the interictal discharges. We analyze whether the estimated spatial modulation of the HRF is a viable biomarker when localizing the ictal onset zone, besides the BOLD spatial signatures themselves.
We use data of twelve patients with refractory focal epilepsy, whom we previously studied in (tousseyn2014sensitivity; tousseyn2014reliable; tousseyn2015correspondence; hunyadi2015prospective). These patients were selected based on the following criteria: 1) they were adults which underwent presurgical evaluation using EEGfMRI, and for which there was concordance of all the available clinical evidence regarding the epileptic focus; 2) subtraction ictal singlephoton emission tomography (SPECT) coregistered to MRI (SISCOM) images were available for all patients, as well as postsurgery MRI scans when patients were seizurefree (international league against epilepsy (ILAE) outcome classification 1–3 (1, completely seizurefree; 2, only auras; 3, one to three seizure days per year auras; 4, four seizure days per year to 50% reduction of baseline seizure days auras; 5, 50% reduction of baseline seizure days to 100% increase of baseline seizure days auras; 6,more than 100% increase of baseline seizure days auras)); 3) interictal spikes were recorded during the EEGfMRI recording session. This study was carried out in accordance with the recommendations of the International Conference on Harmonization guidelines on Good Clinical Practice with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki, for their data to be used in this study, but not to be made publicly available. The protocol was approved by the Medical Ethics Committee of the University Hospitals KU Leuven. For the complete data on the patients’ etiology, and the number of observed IEDs, we refer to Table 1.
patient  gender 

etiology  surgery 




p01  F  L temporal  HS 

3  5  
p02  F  L parietal  FCD 

4  5  
p03  F 

SturgeWeber  
p04  M  R temporal  unknown  
p05  F 

HS 

1  8  
p06  F  R frontal  FCD 

5  2  
p07  F 

DNET 

1  4  
p08  M 

unknown 


p09  F  L occipital  FCD 


p10  F  R temporal  HS  refused  
p11  M 

HS 

1  6  
p12  F  R temporal  CNS infection  refused 
Functional MRI data were acquired on one of two 3T MR scanners (Achieva TX with a 32channel head coil and Intera Achieva with an eightchannel head coil, Philips Medical Systems, Best, The Netherlands) with an echo time (TE) of 33 ms, a repetition time (TR) of either 2.2 or 2.5 s, and a voxel size of 2.6 3 2.6 mm^{3}
. EEG data were recorded using MRcompatible caps with 30 to 64 electrodes, sampled at 5 kHz. The EEG signals were bandpass filtered offline between 150 Hz, gradient artifacts were removed and pulse artifacts were subtracted. The signal of every channel is divided by its standard deviation. Two neurologists subsequently inspected and annotated the EEG signals for IEDs.
The fMRI images were realigned, slicetime corrected and normalized to MNI space, resampled to a voxel size of 2 2 2 mm^{3}, and smoothed using a Gaussian kernel of 6 mm full width at half maximum (FWHM). These processing steps were carried out using SPM8 (Functional Imaging Laboratory, Wellcome Center for Human Neuroimaging, University College London, UK) (friston1994statistical). We refer the reader to (tousseyn2014sensitivity) for a detailed description of these preprocessing steps.
We regress out covariates of no interest from the fMRI data. These include the six motioncorrection parameters, and the average time series in the white matter and the lateral ventricles (cerebrospinal fluid). If necessary, also boxcar regressors are added at moments of substantial scantoscan head movement (larger than 1 mm based on the translation parameters). To reduce remaining nuisance effects, we regress out the first five principal components of the BOLD time series within the cerebrospinal fluid and white matter regions
(behzadi2007component).Subsequently, the BOLD time series are bandpass filtered between 0.008–0.20 Hz using the CONN toolbox (whitfield2012conn).
The dimensionality of the fMRI data is reduced by means of an anatomical parcellation of the brain. The initial images are segmented into regionsofinterest (ROIs) according to the Brainnetome atlas, which consists of 246 parcels in the grey matter(fan2016human). For every ROI, one BOLD time series is constructed as the average of the time series of all voxels within the ROI.
If multiple acquisition runs (within the same recording session) had been done, the EEG and fMRI data of the different runs are temporally concatenated.
In previous work (vaneyndhoven2019semi), we have shown that preenhancing the EEG signals using a datadriven, spatiotemporal filter that is tuned to maximize the signaltonoise ratio (SNR) of IEDs with respect to the background EEG and artifacts, leads to a BOLD predictor that is more performant than many other predictors, including those based on simple stick functions (lemieux2001event; salek2006hemodynamic), ICA (formaggio2011integrating; abreu2016objective) or EEG synchronization (abreu2018eeg). This preenhancement strategy based on multichannel Wiener filters (MWF) has errorcorrecting capabilities and produces an IED representation that improves the localization sensitivity of EEGcorrelated fMRI (vaneyndhoven2019semi).
In brief, the MWF is estimated by first performing timedelay embedding of the multichannel EEG signals , leading to an extended multichannel, multilag signal
(1) 
and subsequently computing the filter coefficients as
(2) 
where is the covariance matrix of the EEG observed during annotated IED segments (), and is the covariance matrix of the EEG outside of IED segments (). For the full derivation, we refer the reader to (somers2018generic; vaneyndhoven2019semi). The EEG signals are then filtered as . Due to the extension with lagged copies of the signals, channelspecific finite impulse response filters are found. Hence, is a set of spatiotemporally filtered output signals, in which IEDlike waveforms are preserved while other waveforms, which are not specific to epilepsy, are supressed^{1}^{1}1Subsampling of the rows of is needed to reverse the timedelay embedding, collapsing it into a multichannel output signal at time only..
We train MWFs for each patient individually, using the same parameter settings as in (vaneyndhoven2019semi): we embed the EEG signals using
positive and negative lags and compute the final filter using the generalized eigenvalue decomposition, which ensures the positive definiteness property of the subtracted covariance matrix in (
2)(somers2018generic).To preserve the intrinsic multiway nature of the data, we represent the preprocessed EEG and fMRI as a tensor and matrix respectively, which are subsequently factorized jointly. This approach differs from the massunivariate treatment in the traditional GLM, where each voxel is treated individually, and only ‘flattened’ EEG time courses can be entered as regressors. Since epilepsy is manifested with considerable variability between patients, we handle the multimodal data of each patient separately.
We adopt a tensorization strategy based on timefrequency transformation of the EEG data to thirdorder spectrograms (time points frequencies channels). After the preenhancement step described in Section 2.3, we create a spectrogram using the Thomson multitaper method (thomson1982spectrum), applied on nonoverlapping EEG segments with a length equal to one repetition time (TR) of the fMRI acquisition. The squared Fourier magnitudes are averaged into 1 Hz bins, from 1 Hz to 40 Hz. Hence, for every EEG channel, we obtain a spectrogram which is synchronized to the fMRI time series. The time points frequencies channels spectrogram, is further normalized to equalize the influence of each channel and each frequency: the mean of each mode1 fiber (holding one time series of squared amplitudes) is subtracted, and afterwards scale normalization is iteratively carried out over the second and third mode to ensure that each channel and each frequency bin contributed the same amount of variance to the data (bro1997parafac).
The average BOLD time series are stacked in a time points ROIs matrix , where ROIs. We normalize each ROI’s time series by subtracting its mean and dividing by its standard deviation.
EEG and fMRI data are acquired simultaneously per subject, and are thus naturally coregistered along the ‘time’ mode. This is captured in a temporal factor matrix that is common between the EEG factorization and the fMRI factorization. However, the electrophysiological changes that are picked up by EEG vary on a much more rapid time scale than the sluggish BOLD fluctuations that (indirectly) correspond to the same neural process. The neurovascular coupling that describes the relation between these two complementary signals can be described by a convolution with an HRF^{2}^{2}2
In this paper, we use the term ‘neurovascular coupling’ to describe the coupling characteristic between EEG and fMRI temporal dynamics, and make the silent assumption that this characteristic is a proxy/surrogate for ‘neurovascular coupling’ as it is understood in neuroscience: the model that describes BOLD changes in response to electrical neural ‘events’, which take the form of local field potentials at the synapses.
.In previous work, we developed a CMTF model in which the HRF itself is parametrically estimated from the data (vaneyndhoven2017flexible), and a matrix multiplication with Toeplitz structure implements the HRF convolution, as proposed in (valdes2009model). In the same paper, we hinted towards an extension based on multiple basis functions to account for the variability of the HRF over brain regions. In the following, we assume that the time course of each EEG source is convolved with an a priori unknown, ROIspecific HRF, which is a superposition of parametrized basis functions, which leads to a modelled contribution of this source to the ROI’s BOLD signal. Hence, in every ROI , the modeled (unscaled) BOLD time course of the th neural source is
(3)  
(4)  
(5)  
(6) 
Here,
is a factor vector holding the time course of the
th EEG source; is an operator that transforms a set of parameters into a full HRF, represented as a vector ; is an operator that transforms into a Toeplitz matrix by populating the main and lower diagonals with the HRF samples (see also A.1); is the weight for the th HRF basis function in the th ROI; is the Toeplitz matrix holding the total HRF in the th ROI^{3}^{3}3The HRF in every ROI does not depend on , and is hence shared between all sources. In (makni2008bayesian; vincent2010spatially; pedregosa2015data), such a constraint has been used to promote robust estimation..This time course is conceptually equivalent to a regressor in the GLM’s design matrix. We treat the HRF parameter sets as unknown variables, which need to be fitted to the data at hand (lindquist2007validity). By parametrizing each basis function, we embed protection against nonsensical HRF shapes, and against overfitting, since the number of parameters to be estimated is greatly reduced compared to the FIR basis in (glover1999deconvolution; aguirre1998variability). We employ a doublegamma HRF, i.e., each HRF basis function is described by five parameters as .
After tensorization, we jointly decompose the EEG tensor and the fMRI matrix into a set of distinct sources.
The thirdorder EEG spectrogram is approximated by a sum of rank1 terms according to the trilinear canonical polyadic decomposition (CPD) (also referred to as Parallel Factor Analysis (PARAFAC)) as in (miwakeichi2004decomposing; marecek2016can; martinez2004concurrent; vaneyndhoven2017flexible). Each rank1 term describes a source (also called ‘component’) in terms of an outer product () of a temporal, spectral, and spatial signature, respectively. Unlike matrix decompositions, the decomposition of a higherorder tensor into a set of sources is unique, up to scaling and permutation ambiguities, without imposing constraints (under mild conditions).
The fMRI matrix is similarly approximated as a sum of rank1 terms. Coupling arises from the temporal signatures , which are shared between the EEG and fMRI factorization. After processing through a hemodynamic system (as described in Section 2.4.3), each source’s BOLD temporal signature is weighted with a spatial signature .
To accommodate additional structured variation in the fMRI data, that is not related to electrophysiological dynamics, we allow a lowrank term to the fMRI factorization which is not coupled with the EEG factorization. We have empirically found that such a lowrank term can capture structured noise, preventing it from biasing the estimation of the parameters which are coupled with the EEG factorization.
The full sCMTF model is then described as:
(7)  
(8)  
(9)  
(10)  
(11)  
(12)  
(13) 
where and are the lowrank approximations; and hold the residuals of both factorizations; describes the CPD model composed of factor matrices , , , which hold the temporal, spectral and EEG spatial signatures in the columns; the HRF matrices are constructed as in (3)–(6); is the fMRI spatial factor matrix; is the HRF basis coefficient matrix; is a rank term to capture fMRIonly structured nuisance; denotes the Hadamard or elementwise product; denotes the Khatri–Rao product (kolda2009tensor).
Note that the coupled part of is described by nonindependent rank1 terms, or equivalently, by rank block terms. Namely, each rank1 term describes the convolution of the th source’s temporal signature with the th basis function, after which a spatial loading with vector is performed; in all ROIs, there is one sourcenonspecific relative weight for the basis function (captured in ), and sourcespecific amplitudes (captured in
). To limit the degrees of freedom, the Khatri–Rao product in (
12)–(13) expresses that the HRF is shared among all sources, which is a constraint that has earlier been used for this purpose (pedregosa2015data). Hence, there are not spatial coefficients, but , i.e., basis function weights and source amplitudes in all ROIs^{4}^{4}4In this way, the Khatri–Rao structure also breaks the curse of dimensionality in the fMRI decomposition if either the number of sources
or the number of basis functions is high (or both).. This model is depicted in Figure 1, omitting to not overload the diagram.We estimate all parameters of the model in (8) and (11) by iteratively minimizing the cost function , composed of two data fitting terms and two regularization terms as in (acar2014structure):
(14) 
(15) 
where the squared Frobenius norm of a tensor is the sum of its squared elements; and denote the Euclidean or norm and the norm or sum of the elements’ absolute values of a vector , respectively; , , and are positive weights; and are vectors which hold the amplitudes with which each source is expressed in the EEG and fMRI data, respectively. The squared Frobenius norms of the residuals promote a good fit of the lowrank approximations to the data, while the regularization terms penalize excessive source amplitudes and promote a parsimonious^{5}^{5}5The sparsitypromoting properties of the LASSO penalty are most useful in the context of an underdetermined system, with more coefficients than observations, e.g. in dictionary learning. Here, the problem is heavily overdetermined, and we do not expect that the amplitudes an go exactly to zero. However, the
penalty is still a useful heuristic to avoid degenerate components in the EEG’s CP decomposition.
model, similar to the groupLASSO method (acar2014structure; yuan2006model). At the same time, the latter penalty also tends to prevent the occurrence of degenerate terms (bro1997parafac). We minimize (14) using the Structured Data Fusion framework in Tensorlab (sorber2015structured; vervliet2016tensorlab), using a quasiNewton method based on a limitedmemory BFGS algorithm, for 50 independent initializations (see A for details regarding the optimization procedure and parameters). After convergence, each set of estimated factors needs to be calibrated to remove certain ambiguities, and model selection must be performed to pick the best solution, with an appropriate (see B for details).We create statistical nonparametric maps (SnPMs) of the obtained spatial signatures to determine which ROIs sources are significantly (de)activated in relation to the found sources (nichols2002nonparametric; waites2005reliable). To this end, we use permutationbased inference, in which the spatial signatures are compared against their empirically derived distributions, which are obtained via resampling of the fMRI data while freezing the other estimated sCMTF factors^{6}^{6}6
Under the null hypothesis of no significant BOLD effect related to the EEG dynamics, the fMRI data may be temporally reshuffled without a significant loss of fit to the EEG dynamics in
.. To account for serial correlations in the fMRI time series, we use the robust waveletbased resampling approach in (bullmore2001colored) to ensure exchangeability and to preserve spatiotemporal correlation structure of the original data in the produced surrogate datasets. For each fMRI dataset and every sCMTF solution, we generate surrogate fMRI datasets using the procedure in (bullmore2001colored). We resample only the adjusted data , i.e., after removing the components which model variation specific to the fMRI data. We perform inference on a pseudo tstatistic, which we compute for every ROI and for every source as follows:construct a local ‘design matrix’ with all estimated temporal signatures as in (3):
find the new ‘betas’ by solving
convert the betas to a tstatistic per source by dividing them by their estimated standard deviation (see (friston1994statistical; poline2012general)).
Through this procedure, we obtain point empirical null distributions for every source and every ROI. We set the significance threshold as to control the familywise error (FWE) rate at , according to the maximum statistic procedure outlined in (nichols2003controlling). That is, for every source , we form the empirical distribution of the maximal tstatistic over all ROIs, and determine sourcespecific thresholds as the 95%percentile (to test for activation) and as the 5%percentile (to test for deactivation). Finally, we obtain statistical maps for all sources by applying these thresholds to the original spatial signatures , which can be considered as the betas of the unshuffled data.
Furthermore, we create a map of the HRF variability over ROIs. For every ROI, we assess how ‘unusual’ the local HRF is, by measuring its calibrated distance in HRF space to all other ROIs’ HRFs. We use two metrics to quantify this (see C for details on the computation).
Extremity is computed as one minus the average of the absolute values of the correlations between a HRF waveform and all other HRFs’ waveforms.
Entropy
of the HRF waveform is computed as the negative logarithm of the conditional probability of the HRF.
Both for the pseudo tmaps as for the HRF extremity and entropy maps, we furthermore limit the inspection to the 20 ROIs with the highest values, if applicable.
An endtoend overview of our pipeline, from data preprocessing up until statistical inference, is depicted in Figure 2.
We use several metrics to quantify the quality of the obtained sCMTF solutions.
We compare the statistical maps with a ground truth delineation of the ictal onset zone (IOZ) to assess their concordance. This ground truth is the manually delineated resection zone for patients that had undergone surgical treatment and that were seizurefree afterwards (vanhoudt2013eeg; grouiller2011or; zijlmans2007eeg; an2013electroencephalography; thornton2010eeg), or otherwise the hypothetical resection zone, based on concordant evidence from multiple modalities other than EEGfMRI (cfr. Section 2.1), for patients that were ineligible for or refused surgery (tousseyn2014sensitivity). The sensitivity for detecting the IOZ is then computed as the fraction of ‘true positive’ cases, which are determined by the presence or absence of significant activation clusters which overlap the IOZ in the spatial signatures . Following the reasoning in (tousseyn2014sensitivity), we do not consider significantly active voxels or regions outside of the delineated IOZ as false positives. Acknowledging epilepsy as a network disorder, such active regions might reflect seizure or IED propagation, despite not being involved in their generation.
Furthermore, we hypothesize that the spatial variation of the HRF over the brain might reveal additional localizing information regarding the IOZ, i.e., based on considerations explained in Section 1, we assume that the HRF in or near the IOZ might be distorted compared to nonepileptic brain regions. We test this hypothesis by assessing whether those regions correspond to high values in the HRF entropy and HRF extremity maps (see 2.6).
Additionally, we inspect the spectral, spatial and temporal EEG signatures of the extracted sources, and we measure whether the spatial fMRI signatures bear any similarity to known networks of restingstate human brain activity (shirer2012decoding).
Table 3 compiles the results of the steps described in B. For each patient, we select the set of sCMTF factors of rank , which best fulfill the criteria. In all cases, we found at least one such a solution, including an IEDrelated component within that solution. Note that sometimes models with different might score well on different (subsets of) the criteria, so the selection of the rank is inevitably ambiguous. In the next section, we analyze the individual set of results for each patient, based on the selected rank, and we analyze the sensitivity of the results to the choice of .
We analyze for each patient the sources which have been estimated via the sCMTF model. Due to space constraints, we discuss the results of two patients in detail in the next subsections, and include complete results for all other patients in the supplementary material. Every time, we show 1) the thresholded pseudo tmaps of the IEDrelated source in the fMRI domain, both for significant activation as for significant deactivation; 2) maps highlighting the ROIs of high HRF entropy and extremity; 3) the temporal profile (timevarying power) , spatial profile (topography) and spectral profile of each source in the EEG domain; 4) the HRF waveforms in the different ROIs, and the HRF basis functions at convergence of the algorithm. We plot maximally 800 s of the temporal signatures, to ensure readability. For ease of comparison, we always overlay the broadband MWF envelope (with an arbitrary vertical offset for visualization only), which is the reference time course for selecting the IEDrelated component (cfr. B.3). For considerations of space, we generally only show the maps of the fMRI spatial signature for the IEDrelated components, but discuss the maps of other components when relevant. We show five axial slices of each map: in each case, we show two slices near the highest and lowest voxels of the IOZ or significant regions of the fMRI spatial signature (whichever lies furthest); if applicable, the middle slice is the crosssection with most overlap between IOZ and spatial signature, and the two remaining slices lie halfway between this slice and the extremal slices; otherwise all three bulk slices are chosen with equal spacing between the extremal slices. We crossvalidate the maps against known resting state networks (RSN) of human brain activity from the Stanford atlas (shirer2012decoding).
We stress again at this point that a subset of the results is prone to errors due to imperfect sign normalization (cfr. B.1). While it is relatively straightforward to unambiguously determine the ‘right’ sign of the EEG signatures, this is more challenging for fMRI. That is, frequently, the polarity of the HRF waveform is ambiguous, and making the ‘wrong’ choice in a voxel (i.e., the HRF which has the opposite effect of the true physical CBF change) immediately leads to the wrong sign of the spatial coefficients in and their pseudo tvalues for all sources . To track the occurrence of this foreseen failure mode, we also investigate the significant deactivations of the sources^{7}^{7}7Alternatively, it is possible to use a pseudo Fstatistic, e.g. the squared pseudo tvalue, to bypass the sign correction altogether. The downside of such an approach is that it is then impossible to distinguish activation and deactivation, which may be meaningful.. Note that we designed the HRF variability metrics so that they are immune to the polarity of the HRFs. Hence, any high score of the HRF metrics can be reliably interpreted. For each case, we separate the twenty waveforms with the highest entropy score, and report how many of those are found in ROIs that overlap with the IOZ, along with the probability (in the form of a pvalue) that this would occur by randomly sampling as many ROIs (under a given fraction of brain that is covered by the IOZ). Hence, this metric is analogous to the positive predictive value (PPV)^{8}^{8}8The PPV is the fraction of positive predictions (in a classification or hypothesis testing framework) which are truly positive, which equals one minus the false discovery rate (FDR)., albeit no rigorous test has been applied in this case.
We analyze the solution with sources, and show the results in Figure 4 and 5. Besides one clear IEDrelated source, there is one other source that is substantially correlated to the reference time course, but with a homogeneous distribution over the head and an unclear spectrum. This may signify that the IEDs do not follow exactly a rank1 structure in the spectrogram, and that they may be nonstationary in time or space (cfr. the argument made for nonstationary seizures in (hunyadi2014block)). The second source’s pseudo tmap had significantly active areas symmetrically in the left and right parietal lobe, much more focalized than the EEG topography. In the EEG time courses, we found indeed IEDlike waveforms at the times of the peaks in the temporal signature. Hence, we suspect that both sources may reflect the onset and propagation of the IEDs to other areas, respectively. Five out of the twenty ROIs with highentropy HRFs overlapped with the IOZ, and a significant finding is that several of them are highly noncausal, i.e., with a positive peak before zero seconds. Figure 5 confirms this, and also shows that the IEDrelated source is significantly active in different ROIs of the IOZ.
We analyze the solution with sources, and show the results in Figure 6 and 7. There is a clear IEDrelated source, and also an artifactual source at Hz, which is also present in other patients. Due to its relatively consistent occurrence, we hypothesize that this artifact is due to the MR acquisition. For example, it may be a remnant of a gradient artifact which is not adequately removed from the data of some channels, cfr. the observation made in (marecek2016can). Surprisingly, this source is significantly active in an extended area in the occipital lobe, overlapping with the visual network. Both HRF metrics correctly reached extreme values at (distinct) ROIs within the IOZ. The pseudot map of the IEDrelated source shows significantly active ROIs that are concordant with the IOZ, and deactivation of a large part of the default mode network. Furthermore, the IEDrelated source’s EEG topography is very consistent with the clinical diagnosis. The fourth source is active in the default mode network, predominantly in the band (cfr. Figure 9). The fifth source had an unclear spectrum, but its temporal signature corresponds to the occurrence of highamplitude IEDs. Its pseudo tmap shows widespread activations over the brain, which did not include the IOZ. We expect that this component captures the propagation of IEDs, after onset near the IOZ, similarly to patient 3.
We provide an overview of the results w.r.t. IOZ detection in Table 2. All results taken together, the sCMTF results allow a correct detection of the IOZ based on the significant IED activation (10/12 cases), significant IED deactivation (6/12 cases), HRF entropy (9/12 cases) and HRF extremity (8/12) cases. All cases are covered by at least one of the metrics, and all patients besides patient 6 had at least two metrics providing correct and complementary localizing info on the IOZ. For nearly all cases, the IEDrelated component’s time course was highly correlated to a reference IED time course, and its spectrum was plausible. In many, but not all cases, this component’s EEG topography was also consistent with the location of the IOZ, though this notion is slightly fuzzy because of the very different spatial domains of EEG and (f)MRI—hence we do not use the term ‘concordant’. Analysis of the spatial, spectral, and temporal signatures, in combination with inspection of the filtered EEG signals, reveals the identity of RSN oscillations and/or artifacts in the majority of cases. For several patients, we found sources that are active in a narrow spectral band near 33–34 Hz. While this likely reflects a technical artifacts as the result of the MR acquisition, we found no concomitant changes at this frequency in the EEG. This may be the result of the normalization procedure which we applied prior to the decomposition: since every frequency bin was given equal importance, even unnoticeable but structured fluctuations at higher frequencies may be captured in a component.
patient 




20 highestentropy ROIs  
ID  activation  deactivation  entropy  extremity  # in IOZ  (pvalue)  
p01  6  no  no  no  yes  yes  1  (0.34)  
p02  3  no  yes  yes  yes  no  1  (0.59)  
p03  2  no  yes  no  yes  yes  5  ()  
p04  4  yes  yes  no  yes  yes  2  (0.32)  
p05  5  yes  yes  yes  yes  yes  6  ()  
p06  2  no  yes  no  no  no  0  /  
p07  4  yes  yes  no  yes  no  1  (0.57)  
p08  2  no  yes  yes  no  no  0  /  
p09  2  yes  yes  yes  no  yes  0  /  
p10  5  yes  yes  no  yes  yes  3  (0.02)  
p11  2  yes  yes  yes  yes  yes  4  (0.01)  
p12  2  no  no  yes  yes  yes  7  () 
For many patients, selecting is ambiguous, since more than one solution (with different ) score well on some of the criteria (cfr. Table 3). Therefore, we analyze impact of the choice of on the sCMTF results. For each patient, we select the solution with the rank which is next in line, i.e., which would be a second best (or equally good) choice, based on the same criteria. This is the solution with for patient 12, for patients 1, 2, 5 and 7, for patients 3, 4, 6, 8, 9 and 10, and for patient 11. For patients 1, 6 and 8, the results deteriorate drastically, as no metric correctly localizes the IOZ. For patient 11, no ROI within the IOZ is significantly activated due to IEDs anymore, but the HRF metrics are still informative. The results for patients 9 and 12 improve, since all metrics are now sensitive to the IOZ. For the other patients, the situation stay more or less the same, i.e., the same metrics are valuable for IOZ localization. However, the maximum value under the different metrics is generally attained at different ROIs compared to the initially selected model.
We have proposed an integrated and structured coupled matrixtensor factorization (sCMTF) framework, which can be used to make inferences on the localization of the ictal onset zone in refractory epilepsy based on simultaneous EEG and fMRI recordings. Our approach aims to perform blind source separation of the neural activity related to interictal epileptic discharges (IEDs), and to characterize it in the spatial, temporal, and spectral domain. To this end, we developed a pipeline consisting of 1) semiautomated EEG enhancement based on annotations of the IEDs; 2) modalityspecific preprocessing and tensorization steps, which lead to a thirdorder EEG spectrogram tensor varying over electrodes, time points, and frequencies, and an fMRI matrix with BOLD time courses for a predefined set of regions of interest or parcels; 3) coupled matrixtensor factorization of the EEG tensor and fMRI matrix along the shared temporal mode, while accounting for variations in the local neurovascular coupling; 4) automated selection of a robust, and relevant IEDrelated component, and nonparametric testing to infer its spatial distribution in the brain.
We have stressed the importance of and accounted for the variability of the hemodynamic response function (HRF) over different patients and brain regions, by equipping the CMTF with the required expressive power via a set of adaptive basis functions. Moreover, after estimating the EEG and fMRI factor signatures, as well as the HRF parameters, we have computed different summary metrics (entropy and extremity) that measure the local deviance of a ROI’s HRF compared to other HRFs in the same brain, and have crossvalidated the spatial map of these metrics against the ground truth localization of the ictal onset zone.
The sCMTF pipeline managed to provide correct detection in all twelve patients in this study, with varying degrees of certainty. The statistical nonparametric map (SnPM) of the spatial signature of the IEDrelated component, obtained with the sCMTF, is the best biomarker The ROIs that are significantly activated under influence of IEDs were the best biomarker in the study, which is in line with the traditional EEGcorrelated fMRI approach (lemieux2001event). In the large majority of patients, several of these regions overlapped with the IOZ. Also the HRF entropy, as a measure of how unlikely an HRF is within a specific set of other HRFs, is a very good biomarker, which almost always identified regions of the IOZ that were complementary to those already found by tracking significant IED activation. In roughly half of all cases, we also found regions within the IOZ that significantly deactivated in association to IEDs. Compared to our earlier work on these data in (vaneyndhoven2019semi), we achieved an additional correct IOZ detection (for patient 11), which is probably thanks to the increased flexibility of the current model. In a waterbed effect, the pseudo tmap for patient 1 no longer allowed correct IOZ detection compared to the simpler pipeline. Luckily, however, this gap is filled here by the new HRF entropy metric. We inspected the 20 HRFs and ROIs with the highest extremity and entropy. Hence, it is inevitable that some or most of these ROIs are not within the IOZ. Standalone HRF metrics would hence have a high false discovery rate, even though for several patients, the high proportion of IOZcovering ROIs among the 20 selected ROIs was very unlikely due to chance (as measured with pvalues). However, the ROIs that were highlighted by the HRF metrics were often distinct from the ROIs identified as significantly activated to the IEDs. Hence, the SnPMs of the IEDs and the entropy metrics provide very complementary information, and when analyzed jointly, they may infer the location of the IOZ with much more certainty, i.e., in the brain area where both IEDrelated and HRFrelated metrics have a high value.
There were substantial differences in (estimated) neurovascular coupling over patients and brain regions, as expected. Since we used ‘regularized’ basis functions, which are parametrized as smooth gamma density functions, the resulting HRFs generally had a plausible shape. However, in some cases we found nonsensical shapes, in which, e.g., the waveform had the same polarity over the whole time course, potentially with a bimodal shape (cfr. patient 4). This serves as a humble reminder to not blindly trust the outputted HRFs (or other factor signatures, for that matter). While we have empirically verified that the optimization algorithm converges properly to the true factor signatures and HRFs for synthetic data under mild conditions, there is no guarantee that this holds true for reallife data, which are orders of magnitude more complex, so that a linear generative model like the sCMTF may not be sufficient to describe the interplay between EEG and fMRI. Moreover, the proper behavior of the sCMTF estimation depends on careful preprocessing, and on a proper selection of hyperparameters (in casu: a good value for the number of sources
). Hence, manual inspection of the data quality and the solution are still required. Even if the estimated HRFs or factor signatures may not fully reflect the ‘correct’ underlying physical phenomena, we have demonstrated that they offer actionable information. Not in the least, via summarizing metrics such as HRF entropy and extremity, our algorithm manages to be reasonably robust to subtle changes in the waveform—which is less of interest here than spatial cues towards the IOZ.The algorithm used its modeling freedom to fit ‘noncausal’ HRFs, which are ahead of the EEG by as much as 10 s. Generally, we indeed found that many of the estimated HRFs had significant positive or negative amplitudes already before the neural correlate visible on the EEG. This is in line with recurrent findings on BOLD changes that precede the IEDs which were observed in the EEG (hawco2007bold; pittau2011changes; jacobs2009hemodynamic; moeller2008changes). We stress that this noncausality may only be in the observation, and not in the underlying physical chain of events: here, it strictly means that we observed BOLD changes in the fMRI data that occur before the corresponding observed neural correlate on the EEG. Despite the fact that many of the HRFs differed substantially from the canonical HRF, which is causal and peaks approximately 6 s after its neural input, we obtained good results as well with the latter HRF as a nonadaptive model for neurovascular coupling (vaneyndhoven2019semi). The reason for this agreement between these different models—which differ substantially in terms of flexibility—is likely that the canonical HRF is positively correlated to the true HRFs which are found inside the IOZ, and as such the resulting activation maps may still be sufficiently informative. In our data and sCMTF results this is indeed the case for many patients.
Importantly, our pipeline heavily relies on a prior enhancement of the interictal spikes in the EEG data, which would otherwise have a too low SNR for the sCMTF algorithm to pick up IEDrelated sources. We employ multichannel Wiener filters, which solely rely on the annotation of a sufficient amount of IEDs in the data itself, or in related data (e.g., data from the same patient, recorded outside the MR scanner). While this task still frequently relies on the skill of human EEG readers and neurologists, advanced automated solutions for interictal spike detection are available (wilson1999spike; scheuer2017spike). Within each solution of a specific rank, we picked the IED component as the one with the highest correlation with a reference time course directly derived from the enhanced EEG. Some of the presented results make clear that this reference time course is not completely free from artifacts, hence caution is warranted when many highamplitude artifacts are still present in the reference. In this study, however, we have not encountered any issues that seemed to be the direct results of a noisy reference during IED component selection.
Overall, the sCMTF pipeline succeeded in extracting meaningful IEDrelated components, alongside components that modeled restingstate neural fluctuations and physiological and technical artifacts. The fact that the sCMTF can estimate signatures and statistical maps for multiple components is a powerful advantage over classical EEGcorrelated analysis. As we demonstrated in the experiments, artifactual influences may be isolated in separate components, which could reduce their impact on IED mapping in the brain. Additionally, we encountered cases where two components were correlated to the IED occurrences: the component with the highest temporal correlation to a reference IED time course then correctly revealed the localization of the IOZ, while the other component presumably modeled the propagation of IEDs to remote brain regions. This observation is analogous to the finding in (hunyadi2016fusion), where a different type of CMTF was applied to average EEG waveforms of IEDs and statistical BOLD maps, which revealed a dissociation between the early IED spike and the subsequent wave, which were related to the onset and spread of the IEDs, respectively. Since we transformed the data with a timefrequency transform that used windows of length TR, our algorithm is unable to unravel different phases within one IED, since they occur in a shorter time frame. However, we identified these different IEDrelated sources by their significantly correlated temporal signatures, and their distinct spatial and spectral profiles. While we did not impose nonnegativity constraints, many estimated EEG spectral and spatial signatures were approximately nonnegative. This need not be the case, however, since the EEG data are normalized in a way such that the resulting signatures would reveal relative increase/decrease, rather than absolute timevarying spectral power (marecek2016can). For example, if a certain component is associated with a power increase in one spectral band, and a simultaneous power decrease in another band, this would be reflected in a spectrum with both a positive and a negative peak.
The endtoend sCMTF pipeline can provide a richer set of results compared to classical EEGcorrelated fMRI analysis. In this respect, it is a more powerful data exploration tool. The tradeoff to be made is that significant computation time goes into the sCMTF and subsequent inference—if one wants to apply it as rigorously as we have done in the current experiments. We seem to be doing a lot of unnecessary work, by computing the sCMTF factors for several numbers of sources, and by repeating the optimization several times for a fixed number of sources. Unfortunately, both ways of repetition seem required to obtain robust results, as we have argued in A and B. However, the EEGonly CP decomposition, which lies at the heart of our initialization strategy, seemed very robust: we found highly similar EEG signatures for almost all random initializations. Probably, this is thanks to the use of the powerful Gauss–Newtontype of optimization. Hence, fewer repetitions of the sCMTF may be already sufficient to arrive at the same robust results. Despite the very reproducible EEG signatures in the initial CP decompositions, we still performed 50 repetitions of the sCTMF, each time slightly varying the initial HRF parameters. As such, we believe our findings are reasonably robust to poor initialization of the HRFs. Performing the sCMTF for many choices of may still be required, as the quality of the result depends on the extraction of an appropriate number of sources. Luckily, however, most ‘optimal’ ranks were quite low in our study, under the heuristic selection procedure. This is reassuring, as it signifies that a typical rank encountered in this context is not problematically high, and no prohibitive computations are needed. Furthermore, we have demonstrated in our experiments that the summary metrics (sensitivity for localizing the IOZ based on different statistical scores) are fairly robust to the choice of , although the estimated signatures themselves differ.
For many patients, the available data was split across multiple runs (i.e., with a few minutes break in between), and we opted to temporally concatenate data over runs, as explained in Section 2.2
. While this violates the coupling model based on HRFs for time samples near the boundaries, we consider the effect minimal, given that the number of those ‘affected’ time samples represents a very tiny fraction of the whole time series. However, a more rigorous approach would be to ‘inverseimpute’ these samples and consider them as missing values: as such, they are ignored during the sCMTF optimization and will not affect the results.
Due to the repeated decomposition and the nature of the nonparametric inference, the computations are highly parallellizable. For a typical dataset with available IED annotations, and with the parameters we have used for this study, the endtoend computation for one patient took a few hours on a machine with twelve cores. To alleviate the computational burden, we have parcellated the fMRI data into 246 regions, based on the Brainnetome atlas (fan2016human)
. This is clearly suboptimal, as the atlas is not patientspecific, and is mostly designed to study healthy brains. There is a serious risk for partial volume effects, in which the IOZ is scattered over several ROIs. As such, the IEDrelated BOLD changes in the part of the IOZ that falls within a certain ROI may get swamped by the remaining BOLD fluctuations within the ROI delineation. Hence, we hope to be able overcome this problem, either by algorithmic improvements, including a speed up of the optimization, or by the use of a patientspecific parcellation or PCAlike compression of the fMRI data. As of yet, it is hard to say whether the fixed atlas had an adverse effect on the results, and it is not so straightforward to compare the statistical maps from this study to maps which are voxelbased. We are currently pursuing experiments in which we employ a hierarchical parcellation: in a first step, the BOLD time series are grouped (but not yet averaged) according to the Brainnetome atlas; subsequently, we use spectral clustering to further refine each Brainnetome parcel based on the correlation matrix of its BOLD time series. As such, this hybrid approach combines a fixed, coarsegrained atlas with a further datadriven subdivision, which can mitigate partial volume effects, while still providing a significant data compression. Alternatively, it is possible to achieve a data reduction while still preserving voxelwise BOLD signals, by limiting the scope of the sCMTF to an a priori defined ROI (e.g., based on a clinical hypothesis stemming from other modalities).
In summary, we have developed and empirically validated a fully integrated framework for EEGfMRI data fusion, which yielded a rich characterization of the interictal activity in time, space, and frequency, and which accounts for and exploits neurovascular coupling variation over the brain. The ability to separate local (de)activation of IEDs from local deviations in the HRF makes the sCMTF a powerful tool for exploratory analysis of interictal EEG and fMRI data. We envision that this approach, with some minor modifications, may also be used to analyze restingstate^{9}^{9}9Since in such a case, no IEDs are present, EEGenhancement like we have done in this study would no longer take place. However, an MWF may still be used to clean up the EEG, e.g. by annotating artifactual periods, which can be removed from the data by the MWF in its dual form (or another tool) (somers2018generic). EEGfMRI activity.
Our complete MATLAB code to execute the pipeline is available at https://github.com/svaneynd/structuredcmtf.
The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/20072013)/ERC Advanced Grant: BIOTENSORS (no. 339804). This paper reflects only the authors’ views and the Union is not liable for any use that may be made of the contained information. This research furthermore received funding from the Flemish Government under the “Onderzoeksprogramma Artificiële Intelligentie (AI) Vlaanderen” programme; from the Bijzonder Onderzoeksfonds KU Leuven (BOF) under the project numbers C24/15/036 and C24/18/097; from the Agentschap Innoveren en Ondernemen (VLAIO) under the project number 150466; from the EU for Horizon 2020 projects 766456, 813120 and 813483; and from EIT for the project SeizeIT (no. 19263).
Declarations of interest: none
In section 2.4.3, we derived the implementation of the convolution with an HRF as a left multiplication of the temporal signatures with a Toeplitz matrix, whose diagonals hold the HRF samples. For a causal convolution, in which the BOLD signal strictly lags its neural correlate, if , otherwise, hence the matrix is lower triangular. This is the situation depicted also in Figure 1.
However, a recurring observation is that BOLD changes can be observed that precede the IEDs themselves (hawco2007bold; pittau2011changes; jacobs2009hemodynamic; moeller2008changes). Hence, we allowed noncausal HRFs that start at most 4 samples before the EEG, which allows for BOLD responses preceding the IEDs by up to 10 s at a typical TR of 2.5 s.
Since the cost function in (14) is nonconvex, any optimization procedure can only guarantee to converge to a local optimum, hence selecting a good starting point is crucial to obtain a reliable solution.
Firstly, we decomposed the EEG data individually according to the CP or PARAFAC model (harshman1970foundations; kruskal1977three; bro1997parafac), to obtain a good initialization for the factors in the sCMTF model. To this end, we used a Gauss–Newton algorithm (cpd_nls with 2000 iterations, 400 conjugate gradient iterations for the step computation, and tolerance on the relative cost function update of , in Tensorlab 3.0 (vervliet2016tensorlab)), which we ran 50 times, from randomly drawn initial factors. We observed that the resulting factors lied very often close together over runs, indicating the algorithm had found a robust solution.
We always employed HRFs, which we manually initialized. To assess whether the eventual sCMTF solution was also robust to the initialization of the HRF basis functions, we used a slightly different set of HRFgenerating parameters in each repetition of the optimization. Figure 10 shows some typical HRF waveforms, which are used to generate the Toeplitz blocks in Figure 1.
From there, we initialized also the fMRI factors in the sCMTF model in (11)–(13). We constructed a flattened ‘design matrix’ in (13) and obtained a rough estimate for as via regression—albeit this does not yet disentangle and . To obtain initializations for the individual spatial factors, we exploit the fact that in every ROI , the Khatri–Rao product of the th columns of and corresponds to a rank1 constraint when folded into a matrix (bousse2018linear; beckmann2005tensorial)
; hence, a rank1 truncated singular value decomposition of the folded
th column of leads to the desired vectors (beckmann2005tensorial), which are further refined via a constrained Gauss–Newton algorithm (bousse2018linear). We approximated the residual of the fMRI data under the initialized (coupled) factors using a rank truncated SVD to capture fMRI nuisances. The parameter was chosen as twice the number of acquisition runs that had been done for a subject.After each of the 50 runs of the initialization procedure, we iteratively optimized (14) with a quasiNewton algorithm (sdf_minf with 1000 iterations, and tolerance on the relative cost function update of , in Tensorlab 3.0 (vervliet2016tensorlab)). Both and were divided by their Frobenius norm, such that afterwards , and we chose so their fit had equal contribution to the cost function. For the regularization penalties, we pick , as in (acar2014structure)
We fitted the sCMTF model to each patient’s data for varying number of sources (rank), i.e., . For each choice of , we ran the optimization procedure 50 times, as explained in A.2. Afterwards, the results need to be aggregated, such that clear conclusions on the sources of interest can be drawn. This involves several steps, which we explain below in chronological order.
Some of the factors that are estimated via minimization of (14) are subject to sign and scale ambiguities, which are inevitable in many BSS contexts (sidiropoulos2017tensor; kolda2009tensor). In the EEG factor model in (8), the factor vectors , and belonging to the same component may be multiplied with arbitrary scaling factors whose product is one, without altering the goodness of fit. Similarly, in the fMRI factor model in (11), sign and scale are exchangeable between corresponding columns of and , and between rows of and rows of . However, to conduct proper statistical inference on each source ’s spatial amplitudes, the elements in must be calibrated, in the sense that it must be possible to compare them across ROIs, and determining the sign is crucial to distinguish local activation from deactivation. Hence, we sequentially fix these ambiguities as follows:
For every component :
and are rescaled to unit norm, and is counterscaled
for both and , the sign is flipped if the sum of squares of the negative elements exceeds that of the positive elements; the sign of is adjusted to preserve the global sign of the EEG rank1 term
if the sign of was flipped, also the sign of is flipped to preserve the global sign of the th fMRI block term
For every ROI :
the local HRF is reconstructed using (6)
the th row of is rescaled and signcorrected as to make unit norm and as to ensure that the HRF’s largest overshoot precedes the largest undershoot; the th row of is counterscaled
To assess the reproducibility of the factors, we use the graphstructured clustering algorithm that we proposed in (vaneyndhoven2019identifying) and briefly summarize here. We represent the factor sets for all 50 repetitions of the fitting procedure as , and use a threshold of 0.85 to construct a binary link matrix that encodes similarities between components from different runs of the optimization (empirically, we found that this threshold led to acceptable cluster definitions in this context). Via lowrank matrix approximation of this link matrix, we then obtained clusters of components that were encountered in varying numbers of repetitions. High cardinality of a cluster is then a sign that the involved component is very reproducible or ‘stable’, since it is part of the factor set upon convergence in many repetitions. We suggest to assign higher trust in such components, as opposed to components in small clusters, which are likely specific to one (potentially poor) local minimum. For the further steps, we condensed each cluster to one of its components, i.e., its centroid. In each cluster, the centroid component is defined as the component which has the largest accumulated similarity with all other components in the cluster^{10}^{10}10By extension, the centroid repetition is defined as the repetition to which the centroid component belongs..
Out of the centroids of the clustered components, we identified the (most) IEDrelated component as the one whose temporal signature was most correlated to a reference time course, which is constructed as the average over channels of the MWF’s output signal’s timevarying (broadband) power. This reference is the BOLD predictor we have proposed for EEGcorrelated fMRI analysis in (vaneyndhoven2019semi), and which provides a good baseline for identifying temporal dynamics that are timelocked to the IEDs.
After the previous steps have been carried out for each setting of , we are left to select the rank whose set of results we proceed with. We heuristically determined an appropriate value by selecting the model which fulfills several criteria:
high core consistency of the EEG decomposition
We compute the core consistency diagnostic (CorConDia) (bro2003new) for the EEG tensor in combination with its estimated factor set from the centroid repetition. The consistency describes how suitable a rank CPD is for the given tensorial data and given factors, and is expressed as a percentage (100% being a very adequate model, and percentages below 70–80% indicating that the model is not appropriate)^{11}^{11}11CorConDia is a popular and robust model selection tool for tensor decompositions (morup2009automatic; acar2007multiway; miwakeichi2004decomposing; liu2016detection; papalexakis2016automatic). To compute it, first the core tensor which is most appropriate (in minimum mean squared error sense) for the given data and CPDderived factors is estimated. Subsequently, CorConDia is computed as the fraction of the core tensor’s sum of squares which is due to offsuperdiagonal elements. When for a given set of factor matrices the CP structure is indeed ideal, the core tensor is superdiagonal and CorConDia equals 100%. Note that for a rank1 model, this notion is meaningless, since the core tensor is a scalar, and CorConDia would trivially be 100% always..
similarity to a reference IED time course
We track the correlation between the IEDrelated component’s temporal signature and the reference temporal signature , as explained in B.3. We expect higher correlations to signify a more suitable model, since generally led to good results in our previous study (vaneyndhoven2019semi).
high significance in the IEDrelated spatial map
We track the highest pseudo tvalue in the SnPM that was created based on the IEDrelated component’s spatial signature (cfr. Section 2.6). A high statistical score indicates a good model fit for the IEDrelated component (abreu2018eeg).
ID  core consistency diagnostic (%) 


maximal tstatistic of 


1  2  3  4  5  6  1  2  3  4  5  6  1  2  3  4  5  6  1  2  3  4  5  6  
p01  100.0  99.6  97.1  90.0  76.8  25  19  14  7  16  20  0.35  0.91  0.38  0.94  0.97  0.97  14.2  20.0  21.5  16.0  7.7  20.1  
p02  100.0  94.9  58.0  19.6  7.2  14  23  27  26  17  18  0.10  0.93  0.95  0.83  0.96  0.88  8.7  17.0  23.4  22.6  17.7  16.9  
p03  100.0  95.6  70.7  29.0  28.4  29  23  24  21  23  19  0.29  0.93  0.92  0.93  0.93  0.93  15.2  16.5  14.9  14.3  13.9  15.0  
p04  100.0  87.0  74.4  37.3  12.3  25  22  15  27  15  20  0.09  0.20  0.61  0.71  0.68  0.81  9.7  6.7  15.4  14.3  11.7  23.4  
p05  100.0  98.7  94.2  93.2  89.7  3  14  14  21  13  17  0.28  0.85  0.00  0.01  0.92  0.01  19.6  5.5  6.5  4.4  9.0  4.6  
p06  100.0  94.6  80.9  33.8  15  21  14  14  11  13  0.55  0.92  0.87  0.80  0.80  0.22  9.9  17.1  9.9  9.6  16.4  36.9  
p07  100.0  97.7  96.7  18  17  12  15  11  12  0.09  0.74  0.80  0.80  0.80  0.22  7.4  10.0  9.4  8.6  11.2  9.9  
p08  100.0  99.7  30.6  12  21  21  11  14  23  0.61  0.95  0.44  0.95  0.94  0.94  13.8  16.0  10.1  20.3  15.5  18.0  
p09  100.0  98.8  95.9  90.1  23.6  21  27  21  15  19  19  0.19  0.66  0.67  0.13  0.48  0.49  14.7  22.5  21.9  10.1  8.0  19.2  
p10  100.0  98.0  95.0  91.7  80.3  23  34  13  24  25  14  0.47  0.29  0.96  0.96  0.96  0.19  11.9  12.9  10.0  8.6  9.9  12.8  
p11  100.0  97.9  78.3  49.1  21  18  12  22  13  12  0.65  0.91  0.50  0.74  0.66  0.67  18.8  12.8  12.0  28.9  18.6  12.6  
p12  100.0  97.3  89.0  59.3  69.4  23  14  14  12  15  12  0.78  0.79  0.60  0.07  0.47  0.50  15.8  15.0  10.9  7.9  10.2  5.5 
The extremity of a specific ROI’s HRF is computed as one minus the average of the absolute values of the Pearson correlation between the HRF waveform and all other ROIs’ HRFs waveforms. I.e., for the th ROI, the extremity is computed as
(16) 
Only the first twenty samples (50 s) are considered. Note that the extremity does not change if the (global, not samplewise) sign (polarity) of one or more HRFs changes.
The entropy of a specific ROI’s HRF is computed as the negative logarithm of the probability of this HRF, conditional on all other ROIs’ HRFs. For example, we first estimated a probability density in HRF space based on all other ROIs’ HRFs, and then evaluated this density at the HRF of the ROI under inspection. From every HRF, we considered the first twenty samples, and then estimated a nonparametric multivariate kernel density in 20dimensional space, by placing a multivariate Gaussian probability kernel at the location of each HRF except one. We made this entropy metric insensitive to the signs of the HRFs, by extending the set of HRF waveforms by their flipped counterparts, and computing the nonparametric density using the resulting HRFs in a leaveoneROIout fashion.
(17) 
in which is a Gaussian kernel distance, which is proportional to
(18) 
in which and are column vectors that store the twenty first samples of the HRFs and , and is a diagonal covariance or bandwidth matrix. We used Silverman’s heuristic to set the kernel bandwidths for each individual dimension, corresponding to one HRF time sample (silverman1986density). I.e., the th bandwidth , which corresponds to the HRF amplitudes at sample , is given by
(19) 
in which is the observed variance (over ROIs) of the HRFs’ amplitudes at the th sample.
Comments
There are no comments yet.