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

EEG-correlated 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, mass-univariate 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 data-driven manner by means of a novel structured matrix-tensor factorization: the single-subject EEG data (represented as a third-order spectrogram tensor) and fMRI data (represented as a spatiotemporal BOLD signal matrix) are jointly decomposed into a superposition of several sources, characterized by space-time-frequency profiles. In the shared temporal mode, Toeplitz-structured 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 EEG-fMRI biomarkers, which can inform the presurgical evaluation of epilepsy. We make all code required to perform the computations available.


Introduction
Refractory epilepsy is a neurological disorder suffered by 30% of approximately 50 million epilepsy patients worldwide (World Health Organization, 2019), in which seizures cannot adequately be controlled by anti-epileptic medication. In the preparation of treatment via resective surgery, interictal epileptic discharges (IEDs) can be localized in the brain with simultaneous EEG-fMRI, which provides a good surrogate for mapping the seizure onset zone (Lemieux et al., 2001;Thornton et al., 2010;van Houdt et al., 2013;Grouiller et al., 2011;Zijlmans et al., 2007;An et al., 2013;Khoo et al., 2017). This mapping is often conducted via EEG-correlated 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 (Gotman, 2008;Lemieux et al., 2001;Zijlmans et al., 2007;Thornton et al., 2010;Salek-Haddadi et al., 2003).
The workhorse for conducting EEG-correlated fMRI analyis has been (Salek-Haddadi et al., 2006)-and will likely continue to be (Poline & Brett, 2012)-the general linear model (GLM) framework (Friston et al., 1994). 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 (Poline & Brett, 2012;Monti, 2011;Lindquist et al., 2009). 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 (Rosa et al., 2010;Murta et al., 2015;Abreu et al., 2018;Van Eyndhoven et al., 2019a). In previous work (Van Eyndhoven et al., 2019a), we addressed this issue by pre-enhancing the EEG signals using a spatiotemporal filter that is tuned to maximize the signal-to-noise ratio (SNR) of IEDs with respect to the background EEG.
We have shown that taking the time-varying 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 (Lemieux et al., 2001;Salek-Haddadi et al., 2006), ICA (Formaggio et al., 2011;Abreu et al., 2016) or EEG synchronization (Abreu et al., 2018). 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 IED-related temporal dynamics to BOLD fluctuations (Friston et al., 1998). However, there is insurmountable evidence that the HRF is not fixed, but varies substantially over subjects (Aguirre et al., 1998), over brain regions (Handwerker et al., 2004), with age (Jacobs et al., 2008), or even with stress level (Elbau et al., 2018). 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 (van Houdt et al., 2013;Bénar et al., 2002;Jacobs et al., 2009;Lemieux et al., 2008;Grouiller et al., 2010). 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 (Lindquist et al., 2009;Lindquist & Wager, 2007;Calhoun et al., 2004;Monti, 2011).
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' (Friston et al., 1998), consisting of the HRF plus its derivative w.r.t. time and its derivative w.r.t. the dispersion parameter (leading to a Taylor-like 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 (Glover, 1999;Aguirre et al., 1998). Other researchers have aimed to find a basis set by computing a low-dimensional subspace of a large set of 'reasonable' HRFs (Woolrich et al., 2004) or by fitting nonlinear functions to given fMRI data (Lindquist & Wager, 2007;Van Eyndhoven et al., 2017). Alternatively, multiple copies of a standard HRF, which differ only in their peak latencies, can be used (Bagshaw et al., 2004). Fi-nally, 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 (Ostwald & Bagshaw, 2011;Caballero-Gaudes et al., 2013). Perhaps surprisingly, the authors of (Caballero- Gaudes et al., 2013) 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 time-invariant 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 (van Houdt et al., 2013). Indeed, in several studies HRFs that deviate from the canonical model were found in regions of the epileptic network (Bénar et al., 2002;Lemieux et al., 2008;Hawco et al., 2007;Pittau et al., 2011;Jacobs et al., 2009;Moeller et al., 2008;van Houdt et al., 2013). Several hypotheses have been postulated to explain this varability, including altered autoregulation due to higher metabolic demand following (inter)ictal events (Schwartz, 2007), vascular reorganization near the epileptogenic region (Rigau et al., 2007), or the existence of pre-spike changes in neuro-electrical activity which are not visible on EEG and which culminate in the IED (Jacobs et al., 2009). 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 (Lindquist et al., 2009;Handwerker et al., 2004;Monti, 2011). Data-driven 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 (Mantini et al., 2007;Marecek et al., 2016). 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 (Valdes-Sosa et al., 2009;Lahat et al., 2015;Calhoun et al., 2009), which is in contrast to EEG-correlated fMRI, where EEG-derived 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 matrix-based techniques like ICA without loss of structure or information (Sidiropoulos et al., 2017;Lahat et al., 2015;Kolda & Bader, 2009;Acar et al., 2007). Tensor-based BSS techniques have been used to mine unimodal EEG data by decomposing third-order spectrograms (channels × time points × wavelet scales) into several 'atoms' (also coined 'components' or 'sources'), each with a distinct spatial, temporal and spectral profile/signature Mørup et al., 2006;Marecek et al., 2016), with successful application in seizure EEG analysis (Acar et al., 2007;De Vos et al., 2007). While a tensor extension of ICA for group fMRI data (in the form of subjects × time points × voxels) exists (Beckmann & Smith, 2005), matrix representations of fMRI remain dominant for single-subject analyses. Coupled BSS techniques can estimate components which are shared between both modalities, providing a characterization in both domains . For example, in (Acar et al., 2017(Acar et al., , 2019Hunyadi et al., 2016;Chatzichristos et al., 2018), multi-subject EEG and fMRI data have been analyzed using coupled matrix-tensor factorization (CMTF), wherein the 'subjects' factor is shared between the EEG trilinear tensor decomposition and the fMRI matrix decomposition. In (Hunyadi et al., 2016), 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 default-mode network (DMN) activity. Also single-subject 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 (Martínez-Montes et al., 2004), 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 subject-specific HRF is co-estimated from the available data, along with the components (Van Eyndhoven et al., 2017).
In this paper, we extend this latter technique in order to account not only for subject-wise 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 data-driven filters from (Van Eyndhoven et al., 2019a), 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.

Patient group
We use data of twelve patients with refractory focal epilepsy, whom we previously studied in (Tousseyn et al., 2014a(Tousseyn et al., ,b, 2015Hunyadi et al., 2015). These patients were selected based on the following criteria: 1) they were adults which underwent presurgical evaluation using EEG-fMRI, and for which there was concordance of all the available clinical evidence regarding the epileptic focus; 2) subtraction ictal single-photon emission tomography (SPECT) coregistered to MRI (SISCOM) images were available for all patients, as well as post-surgery MRI scans when patients were seizure-free (international league against epilepsy (ILAE) outcome classification 13 (1, completely seizure-free; 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 EEG-fMRI 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.

Data acquisition and preprocessing
Functional MRI data were acquired on one of two 3T MR scanners (Achieva TX with a 32-channel head coil and Intera Achieva with an eight-channel 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 MR-compatible caps with 30 to 64 electrodes, sampled at 5 kHz. The EEG signals were band-pass filtered offline between 1-50 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, slice-time 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) (Friston et al., 1994). We refer the reader to (Tousseyn et al., 2014a) for a detailed description of these preprocessing steps. We regress out covariates of no interest from the fMRI data. These include the six motion-correction 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 scanto-scan 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 (Behzadi et al., 2007).
The dimensionality of the fMRI data is reduced by means of an anatomical parcellation of the brain. The initial 79 × 95 × 68 images are segmented into regions-ofinterest (ROIs) according to the Brainnetome atlas, which consists of 246 parcels in the grey matter (Fan et al., 2016). For every ROI, one BOLD time series is constructed as the average of the time series of all voxels within the ROI.
Further customized EEG preprocessing steps are treated in Sections 2.3 and 2.4.
If multiple acquisition runs (within the same recording session) had been done, the EEG and fMRI data of the different runs are temporally concatenated.

Multi-channel Wiener filtering for spatio-temporal EEG enhancement
In previous work (Van Eyndhoven et al., 2019a), we have shown that pre-enhancing the EEG signals using a data-driven, spatiotemporal filter that is tuned to maximize the signal-to-noise 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 (Lemieux et al., 2001;Salek-Haddadi et al., 2006), ICA (Formaggio et al., 2011;Abreu et al., 2016) or EEG synchronization (Abreu et al., 2018). This pre-enhancement strategy based on multi-channel Wiener filters (MWF) has errorcorrecting capabilities and produces an IED representation that improves the localization sensitivity of EEGcorrelated fMRI (Van Eyndhoven et al., 2019a).
In brief, the MWF is estimated by first performing timedelay embedding of the multi-channel EEG signals x[t], leading to an extended multi-channel, multi-lag signal . . .
and subsequently computing the filter coefficients aŝ where R xx = E xx T H = 1 is the covariance matrix of the EEG observed during annotated IED segments (H = 1), and R nn = E xx T H = 0 is the covariance matrix of the EEG outside of IED segments (H = 0). For the full derivation, we refer the reader to (Somers et al., 2018;Van Eyndhoven et al., 2019a). The EEG signals are then filtered asŴ Tx . Due to the extension with lagged copies of the signals, channel-specific finite impulse response filters are found. Hence,Ŵ Tx is a set of spatiotemporally filtered output signals, in which IED-like waveforms are preserved while other waveforms, which are not specific to epilepsy, are supressed 1 . We train MWFs for each patient individually, using the same parameter settings as in (Van Eyndhoven et al., 2019a): we embed the EEG signals using τ = 4 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) (Somers et al., 2018).

Higher-order data representation
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 mass-univariate 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.
2.4.1. Spatio-temporal-spectral tensor representation of EEG We adopt a tensorization strategy based on timefrequency transformation of the EEG data to third-order spectrograms (time points × frequencies × channels). After the pre-enhancement step described in Section 2.3, we create a spectrogram using the Thomson multitaper method (Thomson, 1982), 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, X ∈ R Is×Ig×Im is further normalized to equalize the influence of each channel and each frequency: the mean of each mode-1 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 (Bro, 1997).

Spatio-temporal matrix representation of fMRI
The average BOLD time series are stacked in a time points × ROIs matrix Y ∈ R Is×Iv , where I v = 246 ROIs. We normalize each ROI's time series by subtracting its mean and dividing by its standard deviation.

Neurovascular coupling in the temporal mode
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 .
In previous work, we developed a CMTF model in which the HRF itself is parametrically estimated from the data (Van Eyndhoven et al., 2017), and a matrix multiplication with Toeplitz structure implements the HRF convolution, as proposed in (Valdes-Sosa et al., 2009). 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, ROI-specific HRF, which is a superposition of K parametrized basis functions, which leads to a modelled contribution of this source to the ROI's BOLD signal. Hence, in every ROI i v , the modeled (unscaled) BOLD time course z (r) iv of the r-th neural source is Here, s r is a factor vector holding the time course of the r-th EEG source; H is an operator that transforms a set of parameters θ k into a full HRF, represented as a vector h k ; T is an operator that transforms h k into a Toeplitz matrix H k by populating the main and lower diagonals Figure 1 Structured coupled matrix-tensor factorization (sCMTF) of EEG and fMRI data can reveal neural sources that are encoded in both modalities, as well as capture the varying neurovascular coupling between the electrophysiological and BOLD changes. ((a)) The EEG signals vary over time points × frequencies × electrodes. The resulting third-order spectrogram tensor X is factorized according to (8) into R rank-1 components, which each consist of a temporal signature s r , a spectral signature g r and a spatial signature m r . ((b)) The fMRI data consist of the average BOLD signal in different brain parcels or regions of interest (ROIs), represented in a time points × ROI matrix Y, and are factorized according to (11). Neurovascular coupling is modeled as a convolution of the EEG temporal dynamics with a ROI-specific hemodynamic response function (HRF), as in (11)-(13). In this example, each local HRF is represented as a linear combination (encoded by coefficients b k ) of K = 3 optimized basis functions, each populating a Toeplitz matrix H k which implements a convolution through matrix multiplication with the temporal signatures s r . Afterwards, each smoothed component r is spatially weighted by a signature v r . This is accomplished by the elementwise product b k * v r of the HRF basis function-specific coefficients b k and the component-specific amplitudes v r .
with the HRF samples (see also Appendix A.1); b k,iv is the weight for the k-th HRF basis function in the i v -th ROI; H iv is the Toeplitz matrix holding the total HRF in the i v -th ROI 3 . This time course z (r) iv is conceptually equivalent to a regressor in the GLM's design matrix. We treat the HRF parameter sets θ k , k = 1 . . . K as unknown variables, which need to be fitted to the data at hand (Lindquist & Wager, 2007). 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 (Glover, 1999;Aguirre et al., 1998). We employ a double-gamma HRF, i.e., each HRF basis function is described by five parameters as

Coupled matrix-tensor factorization of EEG and fMRI
After tensorization, we jointly decompose the EEG tensor X and the fMRI matrix Y into a set of distinct sources.
The third-order EEG spectrogram is approximated by a sum of R rank-1 terms according to the trilinear canonical polyadic decomposition (CPD) (also referred to as Parallel Factor Analysis (PARAFAC)) as in Marecek et al., 2016;Martínez-Montes et al., 2004;Van Eyndhoven et al., 2017). Each rank-1 term s r • g r • m r 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 higher-order 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 rank-1 terms. Coupling arises from the temporal signatures s r , 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 v r .
To accommodate additional structured variation in the fMRI data, that is not related to electrophysiological dynamics, we allow a low-rank term to the fMRI factorization which is not coupled with the EEG factorization. We have empirically found that such a low-rank term can capture structured noise, preventing it from biasing the estimation of the parameters which are coupled with the EEG factorization.
3 The HRF in every ROI does not depend on r, and is hence shared between all sources. In (Makni et al., 2008;Vincent et al., 2010;Pedregosa et al., 2015), such a constraint has been used to promote robust estimation. The full sCMTF model is then described as: whereX andŶ are the low-rank approximations; E x and E y hold the residuals of both factorizations; S, G, M describes the CPD model composed of factor matrices S ∈ R Is×R , G ∈ R Ig×R , M ∈ R Im×R , which hold the temporal, spectral and EEG spatial signatures in the columns; the HRF matrices H k are constructed as in (3)-(6); V ∈ R Iv×R is the fMRI spatial factor matrix; B ∈ R Iv×K is the HRF basis coefficient matrix; N, P is a rank-Q term to capture fMRI-only structured nuisance; * denotes the Hadamard or elementwise product; denotes the Khatri-Rao product (Kolda & Bader, 2009).
Note that the coupled part of Y is described by RK nonindependent rank-1 terms, or equivalently, by K rank-R block terms. Namely, each rank-1 term (H k s r ) • (b k * v r ) describes the convolution of the r-th source's temporal signature with the k-th basis function, after which a spatial loading with vector (b k * v r ) is performed; in all ROIs, there is one source-nonspecific relative weight for the basis function (captured in b k ), and source-specific amplitudes (captured in v r ). 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 (Pedregosa et al., 2015). Hence, there are not RKI v spatial coefficients, but (R + K)I v , i.e., K basis function weights and R source amplitudes in all I v ROIs 4 . This model is depicted in Figure 1, omitting N, P to not overload the diagram.
We estimate all parameters of the model in (8) and (11) by iteratively minimizing the cost function J, composed of two data fitting terms and two regularization terms as in (Acar et al., 2014): where the squared Frobenius norm A 2 F of a tensor A is the sum of its squared elements; a 2 and a 1 denote the Euclidean or 2 -norm and the 1 -norm or sum of the elements' absolute values of a vector A, respectively; β x , β y , γ x and γ y are positive weights; λ x and λ y 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 low-rank approximations to the data, while the 1regularization terms penalize excessive source amplitudes and promote a parsimonious 5 model, similar to the group-LASSO method (Acar et al., 2014;Yuan & Lin, 2006). At the same time, the latter penalty also tends to prevent the occurrence of degenerate terms (Bro, 1997). We minimize (14) using the Structured Data Fusion framework in Tensorlab (Sorber et al., 2015;Vervliet et al., 2016), using a quasi-Newton method based on a limited-memory BFGS algorithm, for 50 independent initializations (see Appendix 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 R (see Appendix B for details).

Statistical inference
We create statistical nonparametric maps (SnPMs) of the obtained spatial signatures v r to determine which ROIs sources are significantly (de)activated in relation to the found sources (Nichols & Holmes, 2002;Waites et al., 2005). To this end, we use permutation-based inference, in which the spatial signatures v r are compared against their empirically derived distributions, which are obtained via resampling of the fMRI data while freezing the other estimated sCMTF factors 6 . To account for serial correlations in the fMRI time series, we use the robust wavelet-based 5 The sparsity-promoting 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 λx an λy go exactly to zero. However, the 1 -penalty is still a useful heuristic to avoid degenerate components in the EEG's CP decomposition.
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 sr.
resampling approach in (Bullmore et al., 2001) 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 L = 250 surrogate fMRIỸ (l) datasets using the procedure in (Bullmore et al., 2001). We resample only the adjusted data Y − NP T , i.e., after removing the components which model variation specific to the fMRI data. We perform inference on a pseudo t-statistic, which we compute for every ROI and for every source as follows: 1. construct a local 'design matrix' with all estimated temporal signatures as in (3): iv , ∀ l , 3. convert the betas to a t-statistic per source by dividing them by their estimated standard deviation (see (Friston et al., 1994;Poline & Brett, 2012)).
Through this procedure, we obtain L-point empirical null distributions for every source and every ROI. We set the significance threshold as to control the familywise error (FWE) rate at α = 0.05, according to the maximum statistic procedure outlined in (Nichols & Hayasaka, 2003). That is, for every source r, we form the empirical distribution of the maximal t-statistic over all I v ROIs, and determine source-specific thresholds T (1−α) as the 95%percentile (to test for activation) and T (r) (α) as the 5%percentile (to test for deactivation). Finally, we obtain statistical maps for all sources r by applying these thresholds to the original spatial signatures v r , 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 Appendix C for details on the computation).
1. Extremity is computed as one minus the average of the absolute values of the correlations between a HRF waveform and all other HRFs' waveforms. 2. Entropy of the HRF waveform is computed as the negative logarithm of the conditional probability of the HRF.
Both for the pseudo t-maps as for the HRF extremity and entropy maps, we furthermore limit the inspection to the 20 ROIs with the highest values, if applicable.
An end-to-end overview of our pipeline, from data preprocessing up until statistical inference, is depicted in Figure 2.

Model performance
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  Figure 2 Interictal EEG and fMRI data can be analyzed via structured coupled matrix-tensor factorizations (sCMTF), which reveals both spatial localization of interictal discharges (spikes), and also localized deviations in neurovascular coupling between electrical and BOLD fluctuations. (a) fMRI and EEG data are first separately preprocessed (yellow block). The fMRI data (top row) are structured as a time points × regions of interest (ROIs) matrix, after BOLD time courses are averaged within predefined or data-driven parcels. The EEG data (bottom row) are structured as a channels × time points × frequencies tensor, after the signals are enhanced via a multi-channel Wiener filter (MWF) which is calibrated based on spike annotations, and subsequently undergo a time-frequency transform. (b) The sCTMF of the EEG and fMRI data (blue block) reveals temporally, spatially and spectrally resolved components, and captures spatially varying hemodynamic response functions (HRFs) (cfr. Figure 1). We show the EEG temporal, spatial and spectral signatures in Figures 4a and 6a, and the HRFs in Figures 4b and 6b, for two selected patients. To initialize the sCMTF factors, first a canonical polyadic decomposition (CPD) of the EEG tensor is computed, from which the remaining fMRI factors are initialized. The full sCMTF model is then computed N times, from these N different initializations, and the stability of the resulting factors over runs is assessed. (c) Statistical images are created for the patient's data and the corresponding sCMTF factors (green block). From the sCMTF factors, the spike-related component is picked as the one with the highest temporal correlation to the filtered EEG signals's broadband power envelope. A statistical nonparametric map (SnPM) of this interictal spike-related component is created, revealing coactivated ROIs in a pseudo-t-map (red). For every ROI, the entropy (and also the extremity) of the HRF is computed by assessing its likelihood under the distribution of all other ROIs' HRFs, and a map of this metric is constructed (blue) to reveal localized HRF abnormalities. Both maps can be used to form a hypothesis on the location of the epileptogenic zone, as we show in in Figures 5 and 7 for the two selected patients. In this paper, we validated our technique on a set of patients for which the outcome is known.
concordance. This ground truth is the manually delineated resection zone for patients that had undergone surgical treatment and that were seizure-free afterwards (van Houdt et al., 2013;Grouiller et al., 2011;Zijlmans et al., 2007;An et al., 2013;Thornton et al., 2010), or otherwise the hypothetical resection zone, based on concordant evidence from multiple modalities other than EEG-fMRI (cfr. Section 2.1), for patients that were ineligible for or refused surgery (Tousseyn et al., 2014a). 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 v r . Following the reasoning in (Tousseyn et al., 2014a), 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 resting-state human brain activity (Shirer et al., 2012). Table B.3 compiles the results of the steps described in Appendix B. For each patient, we select the set of sCMTF factors of rankR, which best fulfill the criteria. In all cases, we found at least one such a solution, including an IED-related component within that solution. Note that sometimes models with different R 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 R.

Patient-specific model selection
We show the goodness of fit of the estimated factors for the EEG tensor X and the fMRI matrix Y in Figure 3. Due to the normalization steps which have been applied to the data (cfr. 2.2), the sCMTF operates in a regime of moderately high relative approximation errors.
3.2. Spatio-temporo-spectral profiles of interictal discharges We analyze for each patient the sources which have been estimated via the sCMTF model. Due to space con- Figure 3 Goodness of fit of each patient's EEG tensor X and fMRI matrix Y, for varying choices of the rank R in the sCMTF. Naturally, the EEG approximation error decreases monotonically for increasing rank (intra-patient). For the fMRI data, the fit already plateaus for very low R. This is due to the presence of additional, uncoupled components n q • p q in the fMRI factorization, which can absorb some of the variance when the number of coupled components is low, but which lose their relevance at higher ranks.
straints, 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 t-maps 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 (time-varying power) s r , spatial profile (topography) m r and spectral profile g r 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 s ref for selecting the IED-related component (cfr. Appendix B.3). For considerations of space, we generally only show the maps of the fMRI spatial signature v r for the IED-related 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 cross-section 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 cross-validate the maps against known resting state networks (RSN) of human brain activity from the Stanford atlas (Shirer et al., 2012).
We stress again at this point that a subset of the results is prone to errors due to imperfect sign normalization (cfr. Appendix B.1). While it is relatively straightfor-ward 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 v (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 r iv and their pseudo t-values for all sources r. To track the occurrence of this foreseen failure mode, we also investigate the significant deactivations of the sources 7 . 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 p-value) 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 , albeit no rigorous test has been applied in this case.

Patient 3
We analyze the solution withR = 2 sources, and show the results in Figure 4 and 5. Besides one clear IED-related 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 rank-1 structure in the spectrogram, and that they may be nonstationary in time or space (cfr. the argument made for nonstationary seizures in (Hunyadi et al., 2014)). The second source's pseudo t-map 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 IED-like 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 high-entropy 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 IED-related source is significantly active in different ROIs of the IOZ.

Patient 10
We analyze the solution withR = 5 sources, and show the results in Figure 6 and 7. There is a clear IED-related source, and also an artifactual source at ±34 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 (Marecek et al., 2016). 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 pseudo-t map of the IED-related source shows significantly active ROIs that are concordant with the IOZ, and deactivation of a large part of the default mode network. Furthermore, the IED-related 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 high-amplitude IEDs. Its pseudo t-map 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.

Summary of all patient's results
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 IED-related 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)MRIhence 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.

Sensitivity to model selection
For many patients, selectingR is ambiguous, since more than one solution (with different R) score well on some of (a) Temporal (sr, top), spatial (mr, middle), and spectral (gr, bottom) profiles of the 2 sources in the EEG domain, and reference IED time course (sref, in grey).

Figure 4 (a)
In the selected solution for patient 3 (R = 2), both sources have a temporal signature that correlated strongly to the reference IED time course. The first source modeled the main onset of IEDs and was low-frequency and topographically focal, while the second source was spatially and spectrally diffuse and captured the propagation of IEDs to remote areas (cfr. Figure 8). (b) Five out of the twenty most deviant HRFs were found inside the ictal onset zone (bold lines, p < 10 −4 ). These HRFs had main peaks before 0 s, i.e., they led to BOLD changes before the corresponding EEG correlate of the IED was seen.

Figure 5
The statistical nonparametric maps of the IED-related component (top two rows) and HRF entropy/extremity maps (bottom two rows) of patient 3 show concordance with the ictal onset zone (IOZ). Especially the regions of significant IED activation were accurate, but also five out of the twenty regions with the most deviant (highest entropy) HRFs were found in the IOZ (cfr. Figure 4b). The ground truth ictal onset zone is highlighted in dark gray with a white contour.   The sCMTF solution with R = 5 sources was selected for patient 10. One source's temporal signature is highly correlated with the reference IED time course and is identified as the IED-related source, which has a characteristic low-frequency behaviour and with a frontotemporal topography, consistent with the IOZ location. The second source, which has very narrowband power around ±33 Hz, likely captured an artifact of the MR acquisition. The fourth source captured α activity in the default mode network (cfr. also Figure 9). (b) Three out of the twenty most deviant HRFs were found inside the ictal onset zone (bold lines, p = 0.02).

Figure 7
The statistical nonparametric maps of the IED-related component (top two rows) and HRF entropy/extremity maps (bottom two rows) of patient 10 show concordance with the ictal onset zone (IOZ). IED occurrences were associated with significantly active regions in and near the IOZ (top row), and at the same time with a deactivation in a part of the default mode network (second row). Three out of the twenty regions with the most deviant (highest entropy) HRFs were found in the IOZ (cfr. Figure 6b). The ground truth ictal onset zone is highlighted in dark gray with a white contour. 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 v IED of the IED-related component; 2) the significantly activated and deactivated ROIs in the fMRI spatial signature v IED ; 3) the ROIs with strongly deviating HRF waveforms, as measured with 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 criteria (cfr. Table B.3). Therefore, we analyze impact of the choice of R 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 R = 1 for patient 12, R = 2 for patients 1, 2, 5 and 7, R = 3 for patients 3, 4, 6, 8, 9 and 10, and R = 4 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.

A novel EEG-fMRI data fusion framework
We have proposed an integrated and structured coupled matrix-tensor 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) modality-specific preprocessing and tensorization steps, which lead to a third-order 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 IED-related 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 cross-validated 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 IED-related 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 EEG-correlated fMRI approach (Lemieux et al., 2001). 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 Figure 8 The second component in patient 3 likely captured the propagation of IEDs from the irritative zone, given its relatively large correlation to the MWF envelope (cfr. Figure 4a). The ground truth ictal onset zone is highlighted in dark gray with a white contour. 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 (Van Eyndhoven et al., 2019a), 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 t-map 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 IOZ-covering ROIs among the 20 selected ROIs was very unlikely due to chance (as measured with p-values). 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 IED-related and HRF-related metrics have a high value.

HRFs vary strongly over subjects and brain regions
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 real-life 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 sourcesR). 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 (Hawco et al., 2007;Pittau et al., 2011;Jacobs et al., 2009;Moeller et al., 2008). 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 (Van Eyndhoven et al., 2019a). 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.

Prior EEG signal enhancement aids analysis
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 IED-related 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 (Wilson et al., 1999;Scheuer et al., 2017). 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 high-amplitude 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.

The interpretation of components
Overall, the sCMTF pipeline succeeded in extracting meaningful IED-related components, alongside components that modeled resting-state 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 EEG-correlated 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 (Hunyadi et al., 2016), 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 IED-related 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 abso-lute time-varying spectral power (Marecek et al., 2016). 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.

Practical considerations
The end-to-end sCMTF pipeline can provide a richer set of results compared to classical EEG-correlated 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 Appendix A and Appendix B. However, the EEG-only 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-Newton-type 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 R 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 R, 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 'inverse-impute' these samples and consider them as missing values: as such, they are ignored during the sCMTF optimization and will not affect the results. Figure 9 The fourth component in patient 10 seemed to pick up activity in the Default Mode Network (DMN), predominantly in the α band (cfr. Figure 6a).

Strategies to alleviate the computational demand
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 end-to-end 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 (Fan et al., 2016). This is clearly suboptimal, as the atlas is not patient-specific, 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 IED-related 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 PCA-like 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 voxel-based. 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, coarse-grained atlas with a further data-driven 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).

Summary
In summary, we have developed and empirically validated a fully integrated framework for EEG-fMRI data fu-sion, 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 resting-state 9 EEG-fMRI activity.
Our complete MATLAB code to execute the pipeline is available at https://github.com/svaneynd/ structured-cmtf.

Acknowledgment
The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007(FP7/ -2013 /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 Artificile 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
Since the cost function J 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 X individually according to the CP or PARAFAC model (Harshman et al., 1970;Kruskal, 1977;Bro, 1997), to obtain a good initialization for the factors S, G, M 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 10 −8 , in Tensorlab 3.0 (Vervliet et al., 2016)), 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 K = 3 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 θ k in each repetition of the optimization. Figure A.10 shows some typical HRF waveforms, which are used to generate the Toeplitz blocks in Figure ??.
From there, we initialized also the fMRI factors in the sCMTF model in (11)-(13). We constructed a flattened 'design matrix' D = H 1 S . . . H K S in (13) and obtained a rough estimate for B T V T as U T = D † Y via regression-albeit this does not yet disentangle B and V. To obtain initializations for the individual spatial factors, we exploit the fact that in every ROI i v , the Khatri-Rao product of the i v -th columns of B T and V T corresponds to a rank-1 constraint when folded into a K ×R matrix (Boussé et al., 2018;Beckmann & Smith, 2005); hence, a rank-1 truncated singular value decomposition of the folded i v -th column of U T leads to the desired vectors (Beckmann & Smith, 2005), which are further refined via a constrained Gauss-Newton algorithm (Boussé et al., 2018). We approximated the residual of the fMRI data under the initialized (coupled) factors using a rank-Q truncated SVD to capture fMRI nuisances. The parameter Q was chosen In each repetition, one early-peaking, one mid-peaking and one late-peaking HRF were sampled from probability distributions of the HRF's generating parameters (θ 1 , θ 2 and θ 3 , respectively) that was created by applying some multiplicative noise to baseline parameters (shown as the bold waveforms). The support extends to four negative samples, which gives the model the freedom to fit noncausal HRFs (relatively to the synchronized EEG).
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 quasi-Newton algorithm (sdf minf with 1000 iterations, and tolerance on the relative cost function update of 10 −8 , in Tensorlab 3.0 (Vervliet et al., 2016)). Both X and Y were divided by their Frobenius norm, such that afterwards X F = Y F = 1, and we chose β x = β y = 1 so their fit had equal contribution to the cost function. For the regularization penalties, we pick γ x = γ y = 10 −3 , as in (Acar et al., 2014) Table B.3 For every patient (ID 1-12), the sCMTF model can be fitted for a varying number of sources or rank R. We select a 'good' value for R post hoc, based on four criteria which are checked intra-patient: 1) the core consistency of the EEG tensor decomposition should be high (> 70%); 2) the IED-related source should be found in sufficiently many ( 10) of the 50 repetitions of the estimation procedure; 3) the correlation of the IED-related source's temporal signature with the reference time course, namely the MWF's broadband envelope, is preferably high; 4) the maximal pseudo t-statistic for the IED-related source's spatial signature is preferably high.

Supplement: sCMTF signatures for all patients
We show the EEG spatial, temporal and spectral signatures, HRF waveforms, and fMRI spatial maps of significant IED (de)activation and HRF variability for all patients (except patients 3 and 10, whose results have been analyzed in the main text).

Patient 1
We analyze the solution withR = 6 sources. Figure 1 shows the EEG signatures and HRF waveforms. One of the sources is highly correlated to the MWF reference (in grey), which was already known from Table B.3. This IED-related source had a typical low-frequency spectrum, which is expected for the typical spike-and-wave interictal discharges. The topography is relatively diffuse, although the highest amplitudes are mostly in the left hemisphere. This is in accordance with the lateralization of ictal onset zone (left temporal lobe, cfr. Table 1). There are some noteworthy observations to be made about some of the other components. The fourth has an unusually sharp spectrum, is mainly localized on two nonadjacent center electrodes, and is sustained for a single period of many seconds Hence, this component likely captured an artifact (of yet unknown origin), although we spotted no largeamplitude changes in the EEG itself. Similarly, the third source is only present at one frontal electrode, and exists in a frequency range above 20 Hz. It might represent a muscle artifact, e.g., due to frowning or twitching of some muscles in the forehead. The HRFs of all ROIs are shown in Figure 1b. Two of the basis functions seem to have converged to a very similar waveform, which is an unfortunate possibility if two initial HRFs are too close to the same local optimum in their respective parameters. This reduces * Corresponding author Email address: simon.vaneyndhoven@kuleuven.be, simon.vaneyndhoven@gmail.com the expressive power of the basis set, which is clearly visible, since many ROIs have a nearly identical HRF. One of the twenty ROIs with the highest-entropy HRF overlapped the IOZ, although clearly this HRF (bold line) is not among the most dissimilar waveforms for this patient. This is also visible in Figure 2: both the HRF entropy and extremity maps show a small overlap with the delineated IOZ. Despite the good correspondence in the EEG domain, no significant (de)activation of the IED-component is found inside the IOZ.

Patient 2
We analyze the solution withR = 3 sources, and show the results in Figure 3 and 4. As for patient 1, we found a source which is strongly correlated to the MWF envelope, and which had a mostly low-frequency behavior characteristic for spikes. The topography is mostly uninformative, and does not clearly correspond to the patient's clinical data. The third source is mostly present at both sides of the head, is very sparsely active in time, and has a highfrequency content: this is most likely an artifact due to the neck muscles. Again, there is one of the highest-entropy HRFs which belongs to a ROI in the IOZ. Now, the waveform is clearly resolved from the other HRFs, through the strong initial dip (before 0 seconds). Such a dip is sometimes observed in HRFs, but its underlying physiological mechanism is not yet fully understood. It is possible that this dip reflects altered vascular autoregulation near the IOZ (cfr. the explanation in the Section 1 of the main text), or a rapid depletion in oxygen due to IED generation (before the IED becomes visible on the EEG). Figure 4 furthermore shows that the IED-related component is significantly active in parts of the IOZ, and deactive in others. As mentioned earlier, this deactivation may or may not be due to errors in sign correction. Interestingly, the ROI with the high alteration in neurovascular coupling is distinct from both the activated and deactivated ROIs.

Patient 4
We analyzed the solution withR = 4 sources, and show the results in Figure 5 and 6. There is one source which is mostly correlated to the reference (but not extremely, see also Table B.3). This source had a right-temporal focus, conform the diagnosis in Table 1. The second source illustrates the phenomenon of an erroneous sign exchange between the spatial and spectral profiles. Also one of the HRFs has a negative polarity, which is a failure of the sign correction procedure (in this case, because there is exceptionally no positive overshoot). However, the HRF variability metrics are still interpretable, and indeed two ROIs among the ones with the highest-entropy HRFs overlap with the IOZ. The IED component is significantly active in a tiny portion of the IOZ (cfr. Figure 6). The second source is significantly active in symmetrical parts of the parietal lobe. Given its ongoing fluctuation over time, we hypothesize that this source captures a resting state network (RSN).

Patient 5
We analyze the solution withR = 5 sources, and show the results in Figure 7 and 8. There is a clear IED-related component, with a very high correlation to the MWF reference, a typical spectrum, and an anterior-temporal focus, which corresponds very well to the patient's diagnosis (cfr. Table 1). The fifth source seems present at only one channel, and has spectral harmonic at ±17 Hz and ±34 Hz. One of these peaks is reminiscent of the fourth component in patient 1. As Figure 8 shows, the HRF entropy and extremity prove to be strong biomarkers for the IOZ in this case, and also the significant IED activation and deactivation allow correct localization. In Figure 7, it is clear that some HRFs may still have the wrong sign, which means that the interpretation of 'active' and 'deactivated' is flipped in those ROIs. Hence, regions of significant deactivation are in fact significantly activated. The fourth source had a significant overlap with the auditory RSN, and its spectrum reveals activity in the β band.

Patient 6
We analyze the solution withR = 2 sources, and show the results in Figure 9 and 10. One source is strongly correlated to the MWF, while the other source is likely an artifact, given its very sparse temporal profile. Both sources coincide at one high-amplitude peak, by which we infer that this is probably an artifactual period in the signal. Indeed, when inspecting the original EEG signals, we found high-frequency muscle artifacts at these times. This source also had no significant activation in its spatial map, which corroborates its non-neuronal origin. The IED-related source had a broader spectrum than most other cases, and an uninformative topography. None of the ROIs with high-entropy HRFs is located in the IOZ. The pseudo t-map provides correct localization of the IOZ, however.

Patient 7
We analyze the solution withR = 4 sources, and show the results in Figure 11 and 12. We found a clear IEDrelated component, with a characteristic spectrum and a topography which is backed up by the patient's diagnosis (left anterior-temporal IOZ). The fourth source has a very similar topography and spectrum to the fifth source in patient 5. One HRF inside the IOZ had a high-entropy, and is distinguishable from the others by its very sluggish waveform, i.e., it is smeared out in time, with no sharp over-or undershoot. Also the pseudo t-map provided an accurate localization of the IOZ. Notably, in this patient, the extremity metric misses the deviating HRF in the IOZ (while the entropy metric picks it up). The second source overlapped with the frontal part of the default mode network (DMN), and is active in the α and low β bands.

Patient 8
We analyze the solution withR = 2 sources, and show the results in Figure 13 and 14. We found two components which had correlated time courses. At the time of the peaks, we found higher-amplitude events in the EEG with dubious origin, hence they may or may not be artifacts. One of both components is more strongly correlated to the MWF, and its activation is concordant with the IOZ. The second component shows high overlap with the sensorimotor network. For this patient, none of the IOZ's ROIs had extreme values of either HRF metric.

Patient 9
We analyze the solution withR = 2 sources, and show the results in Figure 15 and 16. In this patient, there is only a moderate correlation of a component with the MWF reference time course. This component's topography (left occipital) agrees with the clinical description, however. The HRF extremity (and not the entropy) is high in a small part of the IOZ. Both the significant IED activation and deactivation allow correct localization as well. The second source seemingly captured high-frequency oscillatory activity in the sensorimotor network, similar to the previous patient.

Patient 11
We analyze the solution withR = 2 sources, and show the results in Figure 17 and 18. The IED-related source had a high correlation with the MWF reference, but an odd bimodal spectrum. Its EEG topography is very consistent with the clinical description. Both HRF extremity and entropy are useful biomarkers for the IOZ. The IED activation and deactivation maps each had a very small overlap with the IOZ. The second source is temporally sparse and captures high-frequency EEG variations, which we identified as muscle artifacts.

Patient 12
We analyze the solution withR = 2 sources, and show the results in Figure 19 and 20. Again we observe an IEDrelated source and a seemingly artifactual source with a spectral peak near 34 Hz. Many of the high-entropy HRFs are highly noncausal, and are associated to ROIs inside the IOZ. Hence, with both HRF metrics, the highest-scoring ROIs provides good localization of the HRF. While there are no significantly active ROIs in the IOZ, there are several significantly deactivated ROIs, which may indicate that the sign standardization was not done flawlessly (cfr. also some of the negative-peaking HRFs for patient 10). Surprisingly, the second source had one significantly active ROI, which overlaps with the IOZ, but which did not match its EEG topography. Hence, the nature of this source remains ambiguous.