Identification of Negative BOLD Responses in Epilepsy Using Windkessel Models

Alongside positive blood oxygenation level–dependent (BOLD) responses associated with interictal epileptic discharges, a variety of negative BOLD responses (NBRs) are typically found in epileptic patients. Previous studies suggest that, in general, up to four mechanisms might underlie the genesis of NBRs in the brain: (i) neuronal disruption of network activity, (ii) altered balance of neurometabolic/vascular couplings, (iii) arterial blood stealing, and (iv) enhanced cortical inhibition. Detecting and classifying these mechanisms from BOLD signals are pivotal for the improvement of the specificity of the electroencephalography–functional magnetic resonance imaging (EEG-fMRI) image modality to identify the seizure-onset zones in refractory local epilepsy. This requires models with physiological interpretation that furnish the understanding of how these mechanisms are fingerprinted by their BOLD responses. Here, we used a Windkessel model with viscoelastic compliance/inductance in combination with dynamic models of both neuronal population activity and tissue/blood O₂ to classify the hemodynamic response functions (HRFs) linked to the above mechanisms in the irritative zones of epileptic patients. First, we evaluated the most relevant imprints on the BOLD response caused by variations of key model parameters. Second, we demonstrated that a general linear model is enough to accurately represent the four different types of NBRs. Third, we tested the ability of a machine learning classifier, built from a simulated ensemble of HRFs, to predict the mechanism underlying the BOLD signal from irritative zones. Cross-validation indicates that these four mechanisms can be classified from realistic fMRI BOLD signals. To demonstrate proof of concept, we applied our methodology to EEG-fMRI data from five epileptic patients undergoing neurosurgery, suggesting the presence of some of these mechanisms. We concluded that a proper identification and interpretation of NBR mechanisms in epilepsy can be performed by combining general linear models and biophysically inspired models.

Introduction

The interest in the concurrent electroencephalography–functional magnetic resonance imaging (EEG-fMRI) method as an important imaging modality in epilepsy surgical planning has increased gradually during the last 20 years (1–3). This method generates whole-brain maps of blood oxygenation level–dependent (BOLD) responses evoked by interictal epileptiform discharges (IEDs), which are then used to locate/identify potential irritative zones (IZs) within the brain. IED-evoked BOLD response analysis is attractive for epilepsy surgery planning, owing to its low invasiveness, accessibility, low cost, and efficiency. The standard clinical protocol for using IED-evoked BOLD signal to demarcate IZs includes (a) concurrent EEG and fMRI recordings while the patient is resting or prompted to sleep; (b) IEDs (spikes and sharp waves) are visually identified by EEG technicians; (c) series of IED onsets are convolved with a canonical hemodynamic response function (HRF, 4) to create a set of regressors; (d) the BOLD signal at each voxel is described as a function of these IED-based regressors via a general linear model (GLM) (4); and (e) statistical parametric maps (t-test and F-test) of the linear coefficients are created to detect positive BOLD responses (PBRs). In advanced clinical protocols, IZ detection with the EEG-fMRI technique is performed with flexible parametric HRF models (5–9) to account for voxel and subject variability. Unfortunately, these parametric HRF models misrepresent atypical BOLD responses frequently observed in certain regions of an epileptic brain, e.g., negative BOLD responses (NBRs), precluding detection of the seizure-onset zones (SOZs) in many patients (10–16). Non-parametric HRF models (17–23) could in principle account for HRF misspecification in epilepsy by sacrificing parsimony, but they are computationally intensive, do not include spatial dependencies, and lack mechanistic foundations. These limitations create opportunities to increase the sensitivity of the EEG-fMRI method in epilepsy.

Initial EEG-fMRI clinical studies have associated SOZs with PBRs as a result of a localized hyperemic response triggered by abnormal neuronal excitability [e.g., the pioneer work by (24)]. More recent data suggest NBRs in SOZs might be caused by local circuit inhibitions during after-spike slow-wave components (13), presumably owing to a profound hyperpolarization of pyramidal cells (25, 26) after a fast spike. In general, this is the most accepted mechanism for the NBR found in many experimental paradigms (27–32). The inhibition can occur in the same active region, but exceeding the excitation, provoking negative changes in the overall neuronal activity, resulting in an NBR. For simplicity, we shall refer to these mechanisms as the enhanced cortical inhibition (ECI). Clinical [see reviews by (33) and (34)] and preclinical (35) studies of focal epilepsy point out to the existence of abnormal decreases in the hyperemic/metabolic ratio during frequent/strong epileptogenic activity, which might be linked to an NBR according to computer simulations presented in this study. We shall refer to this mechanism as altered neurometabolic/vascular couplings (ANCs).

However, not all NBRs might be considered as SOZ candidates during surgical planning. For example, deactivations (or disruption) of normal resting state networks (RSNs), such as the default mode network (DMN), have also been linked to NBRs, a phenomenon reported in epileptic patients during IEDs (11, 36, 37). Here, we refer to this mechanism as neuronal disruption of network activity (NDA). Also, an initial work by Harel et al. (38) suggested that arterial blood stealing (ABS) could cause an NBR in healthy brain areas as a result of decreases in cerebral blood flow (CBF) and volume (CBV) in a region in close proximity to an IED-evoked PBR. A recent computational model (39) demonstrated that ABS is physically possible in the brain vasculature. More recent studies have corroborated experimentally the existence of ABS (40, 41). Therefore, areas with NDA and ABS types of NBRs are not IZs; hence, they should not be considered during the surgical workup.

Classifying these different types of NBRs from the noise fMRI signal must be challenging. However, it is reasonable to expect they have different HRF waveforms, which could be used as fingerprints of the underlying mechanism. To verify and take advantage of these differences for the identification of the mechanisms, it is advised to have biophysical models. Incorporating biophysical model–based discrimination of disparate NBR types in refractory focal epilepsy may significantly improve the accuracy of the EEG-fMRI method to localize/delineate SOZs, thereby increasing success rates of ablative surgery. Windkessel (balloon) models (42) have been utilized in the last decades for statistical inference of BOLD signals (43, 44) due to their parsimonious capabilities to capture most of the features of the HRFs reported experimentally.

In this article, we propose a comprehensive Windkessel-based model to account for these four possible mechanisms underlying NBRs in patients with focal epilepsy. Using the model, we predict a specific HRF waveform for each of the four NBR mechanisms aforementioned. We also investigate if these HRFs are classifiable from noisy fMRI data. To that end, HRFs were fitted using the near-neighborhood exogenous autoregressive (NN-ARx) (45) model. HRF dimensionality was reduced using the principal component analysis (PCA). We subsequently build a machine learning ensemble classifier that uses the first three principal components as features and their corresponding mechanisms as classes. We evaluate the performance of the classifier in predicting mechanisms from their BOLD signals. Finally, we used this method to evaluate the presence of different NBR types in cases of drug-resistant focal epilepsy.

Materials and Methods

EEG-fMRI Data

This is a prospective study duly approved by the Western Institutional Review Board (WIRB #20160218). Parents or legal guardians of 10 patients (9–18 y/o) recruited at Nicklaus Children's Hospital signed a written approved–informed consent. All patients were refractory to pharmacology treatment and exhibited frequent IED. In this context, “frequent” was defined as at least 1 IED per minute. Patients needing sedation or vascular malformations were excluded. In this study, patients exhibiting significant NBRs were only included (n = 5). Demographics and clinically relevant findings are summarized in Table 1. The results of EEG source localization, positron emission tomography (PET), ictal single-photon emission computed tomography (SPECT), subdural implantation (ECoG/sEEG), MRI diagnosis, and pathology were annotated. The MR-compatible EEG system used in this study is not Food and Drug Administration–approved. Therefore, results obtained from the EEG-fMRI analysis were not used in any form during the surgical evaluation of the patients.

TABLE 1

Table 1. Patient clinical information.

We acquired four 10-min trials of simultaneous EEG-fMRI data from each patient. Using the fMRI trials for which the IEDs were better identified from the EEG, we fitted GLM and NN-ARx to estimate IED-related PBRs and NBRs. The IEDs were visually detected and classified into several subtypes by two experts based on their morphology, the semiology of the patient, and other neuroimaging modalities. The different subtypes of IEDs and other types of events were used as different types of inputs (conditions) in GLM and NN-ARx. In addition, motion correction parameters were included as nuisance regressors.

MRI data were collected in a Philips 1.5-T scanner with a 16-channel SENSE Rx coil. fMRI was acquired with a GE-EPI sequence. Each fMRI scan consists of 21 interleaved slices 6 mm thick with a 2-mm gap, in-plane voxel size of 3 × 3 mm, and field of view (FOV) = 204 mm. Flip angle (FA) was 90, repetition time (T_R) = 2 s, and echo time (T_E) = 45 ms. For the purpose of anatomical reference, a high-resolution T1-weighted image was acquired using a spoiled three-dimensional (3D) gradient echo sequence with T_R = 9.7 ms, T_E = 4 ms, and FA = 12. The structural MRIs have 90 to 100 slices, covering the whole brain. In some cases, a T2-weighted 3D image was also acquired with parameters: T_R = 25 ms and T_E = 3.732 ms, FA = 30, FOV = 240 mm, and 160 2-mm-thick axial slices. 3D fluid-attenuated inversion recovery volumes were acquired in sagittal plane using the following parameters: T_R = 4,800 ms, T_E = shortest; FA = 40°; 230 sagittal cuts, matrix: 212 × 185 × 230 mm; FOV = 250; matrix reconstruction isovoxel 0.98 mm. The fMRI volumes were preprocessed using statistical parametric mapping (SPM) (http://www.fil.ion.ucl.ac.uk/spm/). They were corrected for motion artifacts and spatially smoothed with an 8-mm Gaussian kernel. Both smoothed and unsmoothed images were used in GLM analysis to detect significant voxels. Although the minimum variance estimator of GLM may be biased due to non-Gaussian noise (46), the latter was necessary to detect near PBR and NBR, as it was suggested by Goense et al. (47) and Harel et al. (38) for detecting the NBR of vascular origin. The GLM analysis and the results were masked to the gray matter using the SPM segmentation obtained from the anatomical image (48).

EEG was recorded using a 10–10 system 32-channel EasyCap (BrainAmp MR, Brain Product GmbH). To record EEG data simultaneously with the fMRI signal, we used MRI-compatible EEG amplifiers (BrainAmp MR, Brain Product GmbH). EEG signals were sampled at 5 kHz and digitized (0.5-μV resolution) (BrainVision Recorder 1.4, Brain Products GmbH). The majority of EEG electrodes had impedances lower than 5 kΩ. The electrocardiogram (ECG) was measured with an ECG electrode attached to the middle of the back of the patients. To synchronize the EEG with fMRI scans, a trigger marking the beginning of the scans was sent to the EEG recording laptop. To ensure the highest temporal precision, the clock of the laptop was synchronized to the 10-MHz clock of the MR console using a Syncbox (Brain Products GmbH). The following EEG preprocessing was performed using BrainVision Analyzer 2 (Brain Products GmbH). To remove the MR-related artifact, the EEG data were first subsampled to 50 kHz using sinc interpolation to virtually increase its resolution and correct the random phase jittering—of no <0.2-ms resolution determined by the 5-kHz sampling rate—that is present in the scan markers. This phase jittering has significant negative effects in the estimation of the MR gradient artifacts as the latter can change as fast as 0.2 ms. Subsequently, we applied a method for removing the MRI gradient artifacts (49), based on the estimation of an average artifact template. The resultant EEG data were bandpass filtered within 0.5–125 Hz. After marking the R-waves using a semiautomatic tool, we applied a method to detect and remove the effects of the balistocardiogram (50). Finally, we applied ICA based on the Infomax method (51, 52) to remove further artifacts.

As IEDs last ~70–200 ms, the input u(t) is modeled as a train of short pulses. To account for the actual time in which the slice containing the region with significant voxels was acquired, the inputs u(t) were transformed according to $u\left(t+n\frac{T_R}{N_z}\right)$, where n is the position of the slice according to the sequence in which they were acquired, and N_z, the number of slices the whole scan.

Ictal SPECT was performed using a Siemens Multispect 3 gamma camera (Hoffman States, IL). Technetium-HMPAO was used as the radiotracer at a dose of 300 μCi/kg with a minimum dose of 3 mCi and a maximum dose of 20 mCi. PET was performed using a GE Discovery-Dimension ST PET/CT system. FDG was injected at a dose of 140 μCi/kg with a minimum of 1 mCi and maximum of 15 mCi. The stereo-EEG was performed utilizing a Natus Neurology, Natus Medical Incorporated, Excel-Tech Ltd (XLTEK) (Ontario, Canada). The electrode type was a Natus Neuro Grass disposable deep cup electrode (silver chloride, AgAgCl).

The Biophysical Model

The proposed model comprises a two-state dynamic causal modeling component (P-DCM) (53) for principal excitatory cells and inhibitory interneurons, extended to having long-range modulatory excitatory inputs also in the inhibitory population. NDA and ECI types of NBRs are explained by adjusting the time constants of the modulatory synaptic connections in the excitatory and inhibitory population of the extended P-DCM model, respectively. Each brain region has a Windkessel component linked to its neuronal activity through an inducing signaling (i.e., the neurovascular coupling). A viscoelastic non-linear delayed compliance was also included (54). To account for blood stealing effects in brain regions sharing a common supply artery (Option 1, Supplementary Table A1), an inductive element was added to connect their respective Windkessel components (39). For regions not sharing a common supplying artery, the two equations for the CBF become independent. For this particular case, a simplified model proposed by Friston et al. (55) is used (Option 2, Supplementary Table A1). An oxygen-to-tissue transport (OTT) component was used to account for the dynamics of the oxygen extraction fraction and the O₂ concentrations in both tissue and blood (i.e., the neurometabolic coupling) (56). ABS and ANC types of NBRs are explained by fitting those parameters in the model controlling the stealing effect size and the vascular/metabolic imbalance, respectively. The differential equations describing the biophysical model are shown in Supplementary Table A1. Supplementary Figure A1 shows the flow diagram for the physiological mechanisms with their respective state variables (Supplementary Table A2). A graphical representation of all model configurations and defining parameters for each of the four mechanisms is shown in Supplementary Figure A2. Values of the parameters for all these particular situations are summarized in Supplementary Table A3. All the specifics related to the particular cases of the model can be found in the Appendix.

Detecting IZs Using the GLM

The fastest and most widely used method to detect significant BOLD responses is the GLM regression (4):

$\begin{array}{l}y\left(t_i\right)={\displaystyle \sum _{k=1}^{N_u}}({\beta }_{hrf}^{\left(k\right)}\left(h\otimes u_k\right)\left(t_i\right)+{\beta }_{der}^{\left(k\right)}\left(h\prime \otimes u_k\right)\left(t_i\right)\\ +{\beta }_{disp}^{\left(k\right)}\left(d\otimes u_k\right)\left(t_i\right))+{\displaystyle \sum _{r=1}^{N_r}}{\beta }_r{x}_r\left(t_i\right)+\eta \left(t_i\right)& (1)\end{array}$

where {t_i}_{i = 1, …, N}, being N the number of scans and N_u the number of types of inputs {u_k(t)}_{k = 1, …, Nu}, {x_r}_{r = 1, …, Nr} the confounding or nuisance regressors (e.g., the motion parameters), and β_r their effect sizes. h, h′, and d are the canonical BOLD HRF (57), its temporal derivative, and its dispersion, respectively. The noise is prewhitened by applying a first-order autoregressive (AR) model to the signal. To demonstrate that the NBR mechanisms are detectable using GLM, we simulated a set of 10-min BOLD signals with a single type of input (N_u = 1) consisting of a random Poisson train of pulses with average frequency of 2.6/min. This input was the same across all models and trials. We then created a simulated set of N = 300 fMRI scans with T_R = 2 s by adding the simulated BOLD signals to the voxels of a real EPI image. The trials from the same type of model were added to neighboring voxels to form spatial clusters. The amplitude of the simulated BOLD signal in each cluster was multiplied by a Gaussian spatial kernel [full width at half maximum (FWHM) = 2.5 mm] with the maximum at the center of the cluster. The resultant set of images was further corrupted with colored noise, according to a spectral density given by 1/f^p, with 0 < p < 1, to account for biological noise inherent to the BOLD signal but not attributable to the temporal filtering of the HRF (58). No nuisance regressor was included in these simulations, i.e., x_r = 0. Following standard fMRI preprocessing procedures, the simulated fMRI scans were spatially smoothed with a Gaussian kernel of 8-mm FWHM. For each voxel, the vectors of coefficient $\beta ={\left[\begin{array}{c}{\beta }_{can}\&{\beta }_{der}\&{\beta }_{disp}\end{array}\right]}^T$ of the GLM were estimated using SPM (4). Using an F contrast, we selected the voxels where the null hypothesis, I₃β = 0 was rejected with p < 0.05, after correcting for multiple comparisons using the family-wise error criterion, i.e., the significant voxels. I₃ is the 3 × 3 identity matrix, and ⊗ denotes temporal convolution.

Estimation of HRFs for the IZs

Although the responses fitted by the GLM, i.e., the first-order Volterra kernel (59), can reconstruct a large variety of HRFs, we used a more adaptable method to accommodate all possible extreme HRF waveforms, the NN-ARx (45), which has been simplified here to extract the HRFs from fMRI time series from the detected IZs. In this particular case, there is no need for a term describing the near-neighbor effect. Under this assumption, the NN-ARx consists of an AR model with an exogenous input s_i, i.e., the IEDs.

$\begin{array}{ll}{y}_i={\mu }_i+{\displaystyle \sum _{j=1}^p}{\phi }_j{y}_{i-j}+{\displaystyle \sum _{j=1}^r}{\theta }_j{s}_{i-j-d}+{\epsilon }_i& (2)\end{array}$

With y_i = y(t_i) defining the BOLD signal at discrete time instants. The term μ_i accounts for any drift in the BOLD signal, which is modeled by a polynomial. ε_i represents a normally distribute noise with zero mean and variance σ to be estimated from the fMRI data. This model is suitable to deal with the colored nature of the fMRI time series and the spatial correlation between neighboring voxels in the images. A simple recursive method to estimate the linear coefficient φ_j and θ_j can be found in Riera et al. (45). By minimizing the Akaike information criterion, NN-ARx estimates the orders and coefficients of the AR (p), the order of a polynomial modeling the drift, and the delay of the HRF onset (d). With these, the HRF can be constructed (45). By applying the NN-ARx to the unsmoothed fMRI scans, we extracted the HRF in the IZs with significant PBR and NBR from the GLM.

Classification of Mechanisms

Once the IZs are detected and their HRF estimated using the NN-ARx method, it is necessary to classify them according to the type of NBR mechanism. For this purpose, we propose to build a classifier based on a machine learning algorithm. We simulated M = 51 trials of 10-min fMRI signals with T_R = 2 s with random inputs consisting of trains of short pulses Poisson-distributed in time. To account for the entire span of HRF waveforms, we model the possible intraindividual and interindividual variability in the parameters by randomly sampling their values, for each trial, from a uniform distribution within the intervals specified in Supplementary Table A3. For all trials, the NN-ARx–based HRFs were normalized by the maximum of their absolute value. The number of time points of the HRFs was $\frac{T}{T_R}=\frac{32}{2}=16$. The PCA was used to reduce the dimension of each HRF to just its three most relevant components (features). Then, a final matrix 5M × 3 of features was created, with number 5 representing the number of classes (i.e., ECI, NDA, ANC, and ABS). The PBR class was included as reference. This matrix of features and the corresponding vector of classes were used to create a multiclass machine learning ensemble classifier based on support vector machine (SVM) (60, 61). This type of algorithm is a very popular and powerful tool for classification and regression in many of the research fields today (62). We tested the ability of the classifier to differentiate among the five different classes trying with different kernel functions to find the optimal classification.

Results

Predicted Responses for Single Impulses

Figure 1 shows the predicted responses for the mechanisms proposed in this article, and their sensibility to the relevant parameters, following a single IED event. A longer inhibitory recovery creates an ECI type of NBR with smaller amplitude (Figure 1A). The recovery time of the network is reflected in the NBR duration of NDA (Figure 1B). For ANC, a disproportionately high neurometabolic to neurovascular coupling ratio yields NBR (Figure 1C). Regarding the vascular phenomena, the NBR only occurs in the presence of blood resistance in the shared vessel. The higher the resistance, the higher the amplitude of the NBR (39) (Figure 1D). As Supplementary Material, we uploaded the model codes to the public. The folder includes a pdf file with the documentation that contains user instructions (http://web.eng.fiu.edu/jrieradi/NBR-Model/).

FIGURE 1

Figure 1. Simulation of the ECI (A), NDA (B), ANC (C), and ABS (D) mechanisms after a single short pulse. Red (blue) corresponds to PBR (NBR). Each column corresponds to variation in a parameter that significantly determines the NBR waveform. The sensibility of the responses to these parameters is illustrated with different curves corresponding to three different values of the parameters covering the ranges in Supplementary Table A3. The light blue arrow indicates how the NBRs change by increasing the value of the parameters. Besides BOLD responses, we also show other candidate observables in MRI: rCBF and rCBV. Increasing the duration of the inhibitory recovery decreases the amplitude of the ECI mechanism. Expectedly, the longer the recovery time constant in the NDA mechanism, the slower the NBR. We also note that higher neurometabolic coupling gain yields higher NBR amplitudes in the ANC mechanism. Besides, the higher the arterial resistance, relative to the arteriole, in the ABS mechanism, respectively, the higher NBR amplitude.

Detection, Estimation, and Classification of NBR Mechanisms

Figure 2 indicates that by using the GLM it is possible to detect voxels exhibiting different simulated NBR types. The mechanism with the least significance is ABS—because its NBR has the smaller amplitude. The ECI mechanism was not included because of the similarity of the HRF with that of the ABS mechanism.

FIGURE 2

Figure 2. Detection, using the GLM, of simulated spatial distribution of PBRs and three NBRs in a realistic sequence of echoplanar fMRI scans. (A) Three-dimensional glass brain showing the regions where the mechanisms were simulated. (B) F contrast (the 3-order identity matrix) and design matrix of the GLM, used to detect significant voxels. (C) Two-dimensional glass brain showing the F statistics. The plots show the responses predicted by the GLM (color curves) and the adjusted data (black dots), for each mechanism. Besides each plot, the estimated values and confidence intervals of the coefficients, β_can, β_der, and β_disp, of the GLM are shown. PBR (red), NDA (green), ANC (cyan), ABS (blue).

Figure 3 shows the performance of the machine learning classifier based on the SVM analysis of ensembles of simulated HRFs as described in section Materials and Methods. We demonstrate the ability of this classifier to predict new mechanisms using a five-fold cross-validation, i.e., leaving five HRFs out for prediction and using the rest as the training set. Six kernel functions were tested (i.e., linear, quadratic, cubic, fine Gaussian, medium Gaussian, and coarse Gaussian). However, we present here only the results obtained with the coarse Gaussian kernels, which provided the best range of accuracy from a minimum 89% to a maximum 93.7%. Trivially, the PBR response is clearly separable from the NBRs. The ANC is distinguishable from the PBRs, even though its HRF can have a significant positive overshoot. NDA and ANC can be distinguished from each other and from the rest of the NBR types, owing to the prolonged recovery of the former and the fast and bipolar shape of the latter. However, the margin of classification and the confidence of prediction of ECI and ABS are the lowest because of their proximity. This means that it might be difficult to distinguish in some cases, at least merely from fMRI signals. In general, mechanisms were incorrectly classified in ~6.3% of the cases, only among ECI and ABS types.

FIGURE 3

Figure 3. NN-ARx estimation and classification of HRFs. (A) Example of a portion of simulated fMRI time series for the ABS mechanism. The input u(t) is represented with a black trace at the bottom graph. (B) Both positive and negative HRFs estimated using NN-ARx from the example in A. (C) Estimated HRFs from all trials. The black continuous curve represents the average across trials, whereas the black dash curve corresponds to a simulation without noise. To geometrically illustrate their separability, the 3D plot depicts the scores of the first three components of the PCA decomposition of the matrix formed by stacking all HRF as row vectors, after normalizing their amplitudes. PBR (red), ECI (yellow), NDA (green), ANC (cyan), and ABS (blue). (D) Results of the five-fold cross-validation of the SVM classifier with a coarse Gaussian kernel function resulted in 6.3% of misclassification. Center figure shows the confusion matrix. TPR, true-positive rate; FNR, false-negative rate; PPV, positive predictive value; FDR, false discovery rate. As expected, the PBR is distinguishable from all the NBRs. NDA and ANC mechanisms were also perfectly classified, whereas ECI and ABS are the ones that are closer to each other.

NBR Mechanisms Associated With Particular IZs in the Epileptic Patients

In this section, we used the HRF classifier, previously trained with data from the biophysical models, to predict the NBR types in some particular IZs of five patients with refractory focal epilepsy. For each case, we applied the following pipeline: (a) GLM-based detection of voxels significantly correlated with the IEDs (SPM); (b) selection of the region-of-interest (ROI) for the IZ of interest; (c) estimation of the NN-ARx HRF, averaged inside spheres within the ROIs; (d) classification of the mechanisms using the machine learning (SVM) classifier; and (e) estimation of key parameters of the biophysical model associated with the identified mechanism. To obtain confidence intervals for the HRFs, we estimated the empirical distribution of the null hypothesis of no significant response using a permutation test. This was done by estimating the NN-ARx HRFs from 5,000 trials with random order of the IEDs. For each time point, the lower and upper confidence values were the 5 and 95 percentiles of these null HRF distributions, respectively. Not all IZs detected for each patient by the EEG-fMRI technique are discussed in this study. Results from these five patients are introduced only as proof of concept. This part of the study does not aim at clinically validating our methodology, but rather at illustrating its value.

Patient 1 is a 14-year-old girl with partial autonomic evolving to tonic seizures and left hemisphere polymicrogyria. We found a PBR–NBR pair surrounding the anterior parietal artery (Figure 4). Slices showing the thresholded F statistics map built from the estimated coefficients of the GLM using SPM overlaid on the T1-weighted image are presented in panel A. The blue crosshair locates the center of the NBR region in the right postcentral gyrus (PG), whereas the red crosshair locates the center of the PBR region in the right superior parietal lobule (SPL), separated by the postcentral sulcus (another PBR in the left SPL is also shown in this panel). The schematic to the bottom illustrates the ABS mechanism—the regions share the final segment of the anterior parietal artery. The dark gray curve in panel B shows the NN-ARx PBR-HRF and its confidence interval, estimated from the real data, and averaged across the voxels satisfying F ≥ 4.5 within a 10-mm radius sphere with origin in red crosshair in panel A. The light gray curve shows the estimated NN-ARx NBR-HRF and its confidence interval and averaged across the voxels satisfying F ≥ 5 within a 10-mm radius sphere with origin in the blue crosshair in panel A. The light red and blue curves show the unnoisy simulated PBR-HRF and NBR-HRF of the fitted ABS model with the estimated value R_A = 0.17. In panel C, the temporal behaviors of the ABS simulated BOLD in the PBR and NBR regions (with the aforementioned estimated parameters) overlap the time series of the real fMRI and the 54 IEDs (input) used in the NN-ARx estimation. Note that, to detect this NBR/PBR pair, the unsmoothed images had to be used, considerably decreasing statistical significance. This is however the strategy used in Goense et al. (47) and Harel et al. (38) to detect close BOLD responses with inverted polarities. According to the predicted HRF type, this NBR should not be classified as an IZ.

FIGURE 4

Figure 4. Patient #1 shows an irritative zone with ABS-type HRF. (A) Slices showing the thresholded F-statistics map built from the estimated coefficients of the GLM overlaid on the T1-weighted image. The blue crosshair locates the center of the NBR region in the Right Post-central Gyrus; whereas the red crosshair locates the center of the PBR region in the Right Superior Parietal Lobule (SPL). The cartoon at the bottom illustrates the ABS mechanism-the regions share the final segment of the Anterior Parietal Artery. (B) The dark gray curve shows the estimated NN-ARx PBR-HRF and its confidence interval, estimated from the real data, and averaged across the voxels within a 10 mm-radius sphere with origin in red crosshair in (A). The light gray curve shows the estimated NN-ARx NBR-HRF and its confidence interval, estimated from the real data, and averaged across the voxels within a 10 mm-radius sphere with origin in the blue crosshair in (A). The light red and blue curves show the unnoisy simulated PBR-HRF, NBR-HRF of the fitted ABS model with the estimated values R_A = 0.17. (C) Temporal behavior of the ABS simulated BOLD in the PBR and NBR regions (with the above-mentioned estimated parameters), the time series of the real fMRI and the input used in the NN-ARx estimation.

Patient 2 is a 13-year-old girl with generalized and uncinate gyral seizures. An NBR with ECI type of HRF was found in this patient (Figure 5). No lesions were present (A). In concordance with the ictal-SPECT (hyperperfusion, B) and the brain source imaging (EEG-BSI) (C), we found an NBR in the right frontal eye field. sEEG data (C) showed better correspondence with a hypometabolism (interictal PET, B) in the right parietal lobe. Ictal EEG points out to a bifrontal spike/slow wave at 3 to 4 Hz. Our classifier linked this particular NBR-HRF to an ECI mechanism (D). A total of 38 IEDs were used to generate the BOLD signal regressors. The right parietal lobe was removed. The patient was free of seizures for 30 weeks. Thermal ablation of the right posterior cingulate gyrus was performed 7 months later obtaining a reduction of 90% of seizure frequency. According to our hypothesis, the right frontal eye field is an IZ with potential to be the SOZ; hence, seizures could resume. An alternative explanation is that the ECI might be reflecting the presence of an inhibitory mechanism linked to the ictal slow-wave component.

FIGURE 5

Figure 5. Patients #2 shows an irritative zone with ECI-type HRF with a high probability to be the seizure onset zone. (A) Two modalities of fMRI anatomical imaging: left, T1-weighted; right, FLAIR. (B) Two modalities of nuclear imaging: left, interictal PET; right, ictal SPECT. (C) Left, Ictal EEG points out to a bifrontal spike/slow wave at 3-4Hz; right, brain source imaging (EEG-BSI) in Right Tempo-Frontal. (D) Composed panel: left top and center, two slices with the thresholded F-statistics map built from the estimated coefficients of the GLM overlaid on the T1-weighted image with blue crosshair locating the center of one NBR region in the Right Frontal Lobe, left bottom-estimated NN-ARx of this particular NBR-HRF was classified as ECI mechanism, right-stereo EEG (sEEG) electrodes placement is shown as reference.

Patient 3 is a 17-year-old boy with partial seizures in the right frontal insula. There are unarguably several regions with NDA type of NBRs linked mainly to the DMN (Figure 6). All nodes of the DMN, as well as the superior frontal gyri, the middle frontal gyri, the left inferior frontal gyrus (IFG), part of the right IFG, the right fronto-opercular region, and the right caudate nucleus, were highly significantly deactivated. In addition, PBR was detected in the right IFG, which could be one of the foci of the IEDs, based on the semiology of the patient and their proximity to the EEG electrodes used to detect the IEDs (i.e., spikes with highest amplitude in electrode F8). Other types of events were also marked and used as regressors in the linear models. Panel A shows the thresholded F statistics map built from the estimated coefficients of the GLM overlaid on the T1-weighted image. Top left axial slice: green crosshair locating the center of one NBR region in the right lateral parietal node of the DMN—coinciding with the maximum value of the F statistics (bottom left sagittal: green crosshair locating another NBR region in the caudate nucleus; right slices: red crosshair locating the center of the PBR region in the right frontal cortex—presumably in the origin of the IEDs). The approximate location of electrode F8 is shown with a green circle in the right axial slice to illustrate the possible relation of the frontal PBR and the IEDs. The right inset shows a short segment of the preprocessed EEG data where 2 IEDs were identified. In panel B, the dark gray curve shows the estimated NN-ARx PBR-HRF and its confidence interval, estimated from the real data around the region marked by the red crosshair in panel A and averaged across voxels satisfying F ≥ 9.5 within a 7-mm radius sphere. The light gray curve shows the estimated NN-ARx NBR-HRF and its confidence interval, estimated from the real data around the DMN node marked by the green crosshair in the top left slice in in panel A and averaged across the voxels satisfying F ≥ 70 within a 10-mm radius sphere. The red and green curves show the unnoisy simulated PBR-HRF and NBR-HRF, respectively, of the fitted model with the estimated value of the recovery time constant: τ_e2 = 3s. For illustration purposes (panel C), we also show the temporal behavior of the simulated BOLD signal in the NBR region (with τ_e2 = 3s), the time series of the real BOLD signal, and the IEDs (input). To account for the actual relative effect of the IEDs, the amplitude of the input pulses was multiplied by the normalized power of the EEG in F8. The patient underwent a right anterior temporal lobectomy and partial hippocampus/anterior–insular resection. The patient is not yet seizure-free.

FIGURE 6

Figure 6. Patients #3 shows an irritative zone with NDA-type HRF in a brain area of the DMN. (A) Slices showing the thresholded F-statistics map built from the estimated coefficients of the GLM overlaid on the T1-weighted image. Top left axial slice: green crosshair locating the center of one NBR region in the Right Lateral Parietal node of the DMN—coinciding with the maximum value of the F-statistics. Bottom left sagittal: green crosshair locating another NBR region in the caudate nucleus. Right slices: Red crosshair locating the center of the PBR region in the right frontal cortex-presumably in the origin of the IEDs. The IEDs were detected using the EEG signal in electrode F8. The approximate location of this electrode is shown with a green circle in the right axial slice to illustrate the possible relation of frontal PBR and the IEDs. The right inset shows a short segment of the preprocessed EEG data where 2 IEDs were identified. (B) The dark gray curve shows the estimated NN-ARx PBR-HRF and its confidence interval, estimated from the real data around the region marked by the red crosshair in (A). The light gray curve shows the estimated NN-ARx NBR-HRF and its confidence interval, estimated from the real data around the DMN node marked by the green crosshair in the top left slice in (A). The red and green curves show the unnoisy simulated PBR-HRF and NBR-HRF, respectively, for a τ_e2 = 3s. (C) For illustration purposes, we also show the temporal behavior of the simulated neuronal activity and BOLD signal in the NBR region [with τ_e2 = 3s], the time series of the real fMRI and the input. To account for the actual relative effect of the IEDs, the amplitude of the input pulses was multiplied by the normalized power of the EEG in F8.

Patient 4 is a 10-year-old boy with focal seizures and leg pedaling. The patient exhibits ANC-type NBR (cyan crosshair), just in the edge of a tumor in the left frontal cortex (Figure 7). Although the null hypothesis in the significant voxels could not be rejected with a probability corrected by multiple comparisons, this probability was set to a very low value (p < 0.0005), and the minimum cluster size of the significant regions was set to five voxels (by decreasing the minimum size of significant voxels, more significant voxels appear in the upper edge of the lesion). Moreover, the HRF was significant according to the permutation test. This HRF corresponded to 50 IEDs identified in electrode F7. Other IEDs, for a total of 81, were also identified and included as regressors in the linear models. Panel A shows an axial slice of the thresholded F statistics map built from the estimated coefficients of the GLM overlying on the T1-weighted image. The IEDs were detected using the EEG signal in electrode F7 (green circle), which was very close to the area with the ANC type of NBR. The right inset shows a short segment of the preprocessed EEG data where 2 of 50 IEDs were identified. The patient suffers from tuberous sclerosis complex (B). The cyan crosshair—the maximum value of the F statistics—locates the center of the NBR region, in the perimeter, and below one of the patient's tumors. The tumor is highlighted with the yellow circle in the axial slice and the red arrow in the coronal slice of the T2-weighted image. In panel C, the gray curve shows the estimated NN-ARx HRF and its confidence interval, estimated from the real data, and averaged across the voxels satisfying F ≥ 6.5 within a 10-mm radius sphere with origin in the crosshair in panel A. The cyan curve shows the unnoisy simulated HRF of the fitted OTT model with the estimated value of the neurometabolic coupling gain: κ = 0.51 s⁻¹. We also show the simulated temporal behavior of g and the simulated BOLD signal in the NBR region (with κ = 0.51 s⁻¹), the time series of the real fMRI, and the IEDs (input). Also, to account for the actual relative effect of the IEDs, the amplitude of the input pulses was multiplied by the normalized power of the EEG in F7. Our EEG-fMRI results predict the SOZ in the periphery of the tumor. The patient has neither been sEEG implanted nor undergone a surgical procedure.

FIGURE 7

Figure 7. Patients #4 shows an irritative zone with ANC-type HRF in a brain area at the edge of a tumor (TSC). (A) Axial slice of the thresholded F-statistics map built from the estimated coefficients of the GLM overlaid on the T1-weighted image. The IEDs were detected using the EEG signal in electrode F7. The approximate location of this electrode is shown with a green circle to illustrate the possible relation between the NBR region and the IEDs. The right inset shows a short segment of the preprocessed EEG data where 2 IEDs were identified. (B) The patient suffers from tuberous sclerosis complex. The cyan crosshair—the maximum value of the F-statistics—locates the center of the NBR region, in the perimeter and below one of the patient's tumors. The tumor is highlighted with the yellow circle in the axial slice and the red arrow in the coronal slice of the T2-weighted image. (C) The gray curve shows the estimated NN-ARx HRF and its confidence interval, estimated from the real data, and averaged across the voxels within a 10 mm-radius sphere with origin in the crosshair in (A). The cyan curve shows the unnoisy simulated HRF of the fitted OTT model with the estimated value of the neuro-metabolic coupling gain: κ = 0.51 s⁻¹. (D) For illustration purposes, we also show the simulated temporal behavior of the state variable g and the simulated BOLD in the NBR region [with κ = 0.51 s⁻¹], the time series of the real fMRI and the input. To account for the actual relative effect of the IEDs, the amplitude of the input pulses was multiplied by the normalized power of the EEG in F7.

Patient 5 is a 9-year-old girl with electric status epilepticus on sleep and left (C3-P3-FZ) ictal spikes–waves discharges. A very significant PBR was found in the premotor cortex with high probability to be the SOZ (Figure 8). There is a non-enhancing cystic lesion in the left posterior frontal lobe (A). Interictal PET (C) reveals hypermetabolism (perhaps due to the high frequency of IEDs) on the left posterior paracentral (both precentral and postcentral sulcus) in agreement with the PBR. EEG-BSI indicates brain sources on the bank of the left central sulcus (C). Hence, we expect total seizure control if this area is resected. In contrast, we found an NBR posterior to the cyst that was classified as ABS. We found a PBR nearby this deactivation, but it was not significant. To illustrate the usefulness of the NN-ARx method, we compared HRFs estimated with it and those obtained with the impulse response function (IRF) method (SPM software). The IRF method was not able to capture underlying HRFs. Results from the EEG-fMRI analysis are shown in panel D.

FIGURE 8

Figure 8. Patients #5 shows an irritative zone with ABS-type HRF. (A) Three different slices of anatomical fMRI imaging (T1-weighted) show non-enhancing cystic lesion in the Left Posterior Frontal Lobe (red crosshairs). (B) Slices showing the thresholded F-statistics map built from the estimated coefficients of the GLM overlaid on the T1-weighted image; left panel, NBR highlighted with blue crosshairs; right panel, PBR highlighted with red crosshairs. (C) Composed panel: top, two slices of nuclear imaging (interictal PET); bottom, EEG-BSI data. (D) Comparison between the NN-ARx method with the impulse response function (IRF) method (SPM software). The IRF method was not able to capture underlying HRFs. Results from the NN-ARx were classified as the ABS mechanism.

Discussion

The most accessible, inexpensive, and least-invasive brain imaging method available for localizing seizure foci with high sensitivity (85%) is fMRI combined with EEG. Unfortunately, its benefits remain controversial for many patients, owing to an incomplete understanding of the neuronal, vascular, and metabolic responses in epileptic tissues, decreasing the sensitivity of the biomarker used to locate foci. A large percentage (~17.4%) of discordance in the use of the EEG-fMRI technique is due to a poor classification of clinically relevant NBR responses (Table 2). Here, we hypothesize that NBRs during IEDs could be caused by both clinically and non-clinically relevant mechanisms. A clinically relevant mechanism should be that resulting as a direct consequence of epileptogenic tissues, i.e., an enhanced inhibition and a vascular/metabolic balance mismatch both due to tissue overexcitability. A secondary effect, such as blood flow stealing and resting-state network shutdown, should be considered not clinically relevant and hence not discussed during the surgical workup. Therefore, tools aiming at the classification of these four mechanisms might increase accuracy in the localization of foci for neurosurgical excision, improving success rates. Henceforth, we discuss the rationale and implications of the mechanisms proposed for NBR genesis in epilepsy.

TABLE 2

Table 2. Results from two studies (different laboratories) about the accuracy of the EEG-fMRI technique, with specifications to (a) typical percentage of refractory epileptic patients who will not be able to complete a successful EEG-fMRI study and (b) typical percentage of these patients with a “discordant” NBR result.

Modeling NBRs With a Neuronal Network Origin

Inhibition-related phenomena, i.e., ECI and NDA, require accounting for the imbalance between the local activation of inhibitory and excitatory neuronal states, which depend on their respective connectivity structure. A two-state model (P-DCM) was used by Havlicek et al. (44, 53) to explain NBR during static and flickering visual stimulation. In this article, we extended the P-DCM model to include an additional external input to the inhibitory population in each brain region of interest and an IED-evoked synaptic modulation of the RSNs. Long-range excitatory (thalamocortical/corticocortical) inputs targeting inhibitory populations in the granular layers of the cerebral cortex have been extensively reported in previous literature. IEDs are mostly initiated by a brief increase in excitatory feedback gains and decrease in the thresholds for firing (63). In many cases, local neuronal excitability is followed by an enhancement in cortical inhibition (e.g., the wave component in the spike-wave events), which has been linked to a robust hyperpolarization in III/V layer pyramidal cells (25, 26). Data by Pittau et al. (13) suggested this type of enhanced inhibition might cause NBR in, or near, the actual IZs. Therefore, we hypothesize NBR with an ECI-HRF type should be included as a potential candidate for ablation in the epilepsy surgical workup. We believe the NBR results from an abnormal enhancement in the external input to inhibitory populations in the neocortex, which is modeled by a large response time of the inhibitory population τ_i1. In this article, the τ_i1 parameter was fitted using the BOLD data to accurately represent the particular NBR waveform. On the other hand, the recovery of the neuronal activity of the disrupted RSN (NDA) after an IED affects one of its nodes is characterized by a modulation of the excitatory synapses in specific areas within the RSNs, which was characterized by a reduced intralaminar excitatory connectivity c_e2 = 0.01 and a large response time τ_e2. The latter actually depends on the way the different nodes interact to effectively “shut down” and recover the network. Note that brain dynamics operate near criticality (64–66), i.e., on the brink to instability. Neuronal activities in this situation require higher recovery time to reach equilibria after perturbed and are associated with large-scale dependencies and scale invariance (67). Therefore, the time response for excitatory τ_e2 was fitted to the BOLD data to accurately characterize the NDA type of NBRs.

Modeling NBRs With a Vascular/Metabolic Origin

Enhancements in the neurovascular coupling gain ε cause increases in CBF, hence a larger PBR effect. Ictal hyperperfusion has been observed with ¹⁵O-H₂O PET and ^99mTc-HMPAO/ECD (33, 34). Using invasive recordings from a preclinical model of epilepsy, we have reported increases in the perfusion gain around the SOZs (35). In this previous study, we associated a small value of κ in the SOZ with an increase in the baseline O₂ metabolism, which might be related to reported glucose hypermetabolism in the SOZs from ictal FDG PET. The interplay between O₂/glucose metabolism and blood perfusion during IEDs is still controversial. Several studies have shown a hypometabolism in IZs, whereas others have reported complex glucose metabolism patterns with hypermetabolism also in some SOZ candidates (34, 68). Using ¹⁵O-H₂O PET, Bittar et al. (68) showed an increase in blood perfusion during IEDs. However, reductions in perfusion have been also reported in the past (34). If this neurovascular coupling gain is kept constant, the parameter that is highly correlated with NBR amplitude will be κ. A disproportionate increase of this value leads to the ANC type of HRF. In this mechanism, the NBR can be seen as an exaggerated initial dip as a result of an abnormally enhanced O₂ metabolism. It is hypothesized that this particular type of NBR may be clinically relevant while defining IZs. NBR can also have a pure vascular origin via an ABS effect. Here, we use a model proposed by Suarez et al. (39) that couples two Windkessels by a common artery to classify ABS types of HRFs. Simulations indicate that the parameter that determines the NBR amplitude is the resistance of the vessel (artery), relative to the total steady state resistance of the vasculature within the tissue, i.e., the arterioles, capillaries, and venules. A vascular anatomical network (VAN) model proposed earlier by Boas et al. (69) predicts also a relative decrease in CBF and O₂ saturation around a brain area undergoing a positive functional hyperemic response. However, the authors are not aware of the application of the VAN model to study NBRs in epilepsy. The existence of ABS effect in the brain have been experimentally demonstrated by several groups (38, 41). Here, we recommend excluding IZs with ABS types of NBRs as potential candidates of SOZ. Some authors have suggested that a type of NBR might also result from the combination of blood backpressure and neuronal inhibition (47, 70). That would explain the presence of the poststimulus overshoot, which is in contrast to those observed in our pure vascular simulations. This undershoot is shown to be determined by the dynamics of the inhibitory neuronal state (44). Other causes including vein delayed compliance might also explain these transient. However, these last mechanisms were not investigated in this study.

NBR Classification

In this article, we focus on the possibility of detecting and classifying the NBR mechanisms using the HRFs extracted from BOLD fMRI signals. We use linear models for the detection and estimation of the HRFs, i.e., GLM and NN-ARx methods, respectively. Under the assumption of the extended balloon model, the validity of the GLM was previously investigated by quantifying the effect size of second-order Volterra kernels (55). Under the same model, the validity of the NN-ARx to estimate HRFs was evaluated in (17, 45). Although linear models have been used to detect and reconstruct BOLD responses for decades, even before addressing the non-linear characteristics of BOLD signals (59, 71), we decided to analyze if linear models are able to accurately characterize PBRs and NBRs in epileptic patients. In general, for balloon/Windkessel models—with sporadic IED events and typical canonical-like responses (57), these linear models are suitable for spatial detection estimation of the HRF. We investigated if the NBR mechanisms can be solely classified from their BOLD responses. This is important for clinical applications when only standard fMRI paradigms are available or designed. Figure 3D shows that the machine learning classifier is able to differentiate PBR (red), ECI (yellow), NDA (green), ANC (cyan), and ABS (blue) in 100% of the cases. However, ECI and ABS were undistinguishable in some trials. In practice, this might be worse as there is a loss of sensitivity and specificity related to the usual misclassifications of IEDs, which is expert dependent. Furthermore, we cannot outline the possibility of having more than one NBR mechanism in the same area at the same time. Thus, failing to identify IZs with multiple NBR mechanisms could lead, in the worst-case scenario, to the incorrect clinical assessment. According to our hypothesis, an ABS/ECI misclassification will be the most critical case. However, these mechanisms are different in nature, and our model predicts different responses when using other imaging or recording modalities. At the expense of experimental feasibility, other imaging technique can be combined with our EEG-fMRI methodology to verify the predicted NBR mechanisms. For example, MION (38) and/or VASO (47) can be used to measure CBV concurrent with BOLD signals. Measurements of neuronal activity can be incorporated using ECoG/sEEG and EEG (27). In addition, CBF can be also included using arterial spin labeling (ASL) (30) and/or FAIR (47). In the case of ABS/ECI misclassification, a radiologist could use EEG concurrently recorded with ASL focused on the particular regions to verify whether there is a decrease in CBF. If no decrease in CBF is observed, ECI mechanisms should be expected. With these multimodal observations, we foresee a significant increase in the margins of classification of the NBR mechanisms, even in the case more than one is present in the same region (47, 72). It is important to highlight that none of the available imaging modalities (Table 1) provides conclusive results, and a thorough data evaluation in the surgical workup is needed for each clinical case. Our approach only aims at providing another level of EEG-fMRI data interpretation to improve the accuracy of this technique. To illustrate this, five specific clinical cases are discussed below.

Epileptic Cases Discussion

We found an ABS type of NBR in patients 1 and 5. In patient 1, an NBR was located in the SPL and a PBR in the PG. The NBR was classified as ECI, which would occur via either U fibers or pure vascular phenomena. However, we cast doubt on ECI as we believe an inhibitory pathway from SPL to PG is rather weak. Note that the PG hosts the primary somatosensory cortex (S1) (73), which is a granular cortex that mainly receives somatotopic feedforward afferents from the ventral posterolateral and posteromedial relay nuclei (VPL and VPM) of the thalamus (74). In addition, the SPL, involved in transforming visual information in complex motor planning, has efferent pathways mainly to the premotor supplementary motor cortices in the precentral gyrus. A top–down inhibition from higher areas (prefrontal cortex) to the somatosensory area is mainly via efferent pathways. Moreover, both difussion spectrum imaging (DSI)-based connectivity (75) and cortical thickness–based connectivity (76) between SPL and PG are rather low. The mechanism might be ABS. Note that the detected BOLD responses are in the vascular domain of the middle central artery, at both sides of the postcentral sulcus. Thus, they might be sharing a final segment of the anterior parietal artery. Although we do not discard the existence of a venous blood backpressure effect, in which the regions could be sharing some anastomotic vein feeding the central sulcal vein or a branch of the superior anastomotic vein of Trolard, it has been reported that venous CBV changes are relevant only for longer stimuli (77). The NBR in patient 5 was found posterior to the cyst. We used data from this patient to illustrate the HRF estimation with the IRF method (SPM) and our NN-ARx method. Because of a probable revascularization around the cyst, angiography data from this patient will be required for a discussion about possible scenarios for the ABS effect. A frontal eye-field NBR with an ECI HRF was found in patient 2, which is most likely due to ictal propagation with frontal slow-wave responses. Slow-wave discharges have been found associated with NBR (13). DMN deactivations, like those found in patient 3, have been systematically reported in the literature for temporal lobe epilepsy (TLE) (11, 37) and for other types of focal epilepsy (36, 78). The pattern of deactivation depicted in Figure 6, with a predominance in the parietal node of the DMN ipsilateral to the focus, is similar to that reported ibyn Fahoum et al. (36) for four of five patients, with concomitant decrease of electrophysiological activity. Our results suggest that the PBR region might have afferents on the anterior part of the caudate nucleus that relay to central nodes of the DMN. This is consistent with the hypothesis of widespread secondary inhibition of non-seizing cortical regions via basal ganglia (79). The temporal profile of the PBR HRF in the right IFG experienced an unpredicted decay (or rebound) correlated with the amplitude of the NBR in the DMN nodes. This might be seen as an interruption of the PBR mechanism by inhibitory afferents coming from the regions exhibiting NDA, which in this case are present all around the right IFG. This might have implications in the interpretations of BOLD responses during IEDs or stimulation paradigms. If the location of the PBR is close to an affected RSN node, its HRF waveform might be misleading of the actual underlying PBR mechanism, due to either the interaction between mechanisms, i.e., inhibitory inputs from the NBR to the PBR region, or the effect of the BOLD spatial point-spread function. This rebound can be also explained as an increase in dHb due to an increase in neuronal activity in the NBR region, or even a vascular reallocation phenomena, as suggested by Hu and Huang (40). They observed positive and negative optical responses, concurrent with local field potentials (LFP) and multiunit activity (MUA) measurements, in rats during hindlimb electrical stimulation. Finally, it has been hypothesized that NDA is a disruption of RSN provoking a reduction of consciousness and cognitive reserve (36). Interestingly, our results suggest that a recovery from this disrupted state is not instantaneous. In our data, the estimated value for the recovery time constant was τ_e2 = 3s. The Epilepsy Connectome Project (ECP) (80) constitutes a huge database that contains clinical, neurophysiological, and resting-state fMRI (rs-fMRI) data of 105 patients with TLE and 55 healthy individuals as control. Using graph (nodes and edges) theory combined with rs-fMRI measures, a characterization of abnormal patterns in the local and global neuronal connectivity in TLE has been possible, thanks to the ECP. Our model-based method to identify different types of NBRs in epilepsy can help provide a neurophysiological foundation to the reported connectivity maps abnormalities. The ANC mechanism, as reported for patient 4, is related to a decrease in the CBF/CMRO₂ balance. NBRs in the hippocampus of rats during bicuculline-induced generalized tonic–clonic seizures were associated with this type of mechanism (81). The authors reported that, even with higher LFP/MUA activity in the hippocampus, as compared to the cortex, the CBF was lower, and the CMRO₂ was higher, yielding to NBRs in the hippocampus. Quantitatively, the unbalance corresponds to a decrease in the ratio $\frac{\epsilon }{\kappa }$ in the OTT model (56), which is $~\frac{0.4}{0.05}=8$ for normal positive responses. Song et al. (35) estimated a ratio approximately five-fold smaller in rat with focal cortical seizures. Although the rCBF/CMRO₂ coupling was reported to be preserved in human IEDs without any apparent lesion (82), we do not discard the possibility of an unbalance produced by a more critical state of tissue pathology. For example, the typical calcification of the surrounding blood vessels present in TSC tumors (83) could hamper the expected IED-induced increase of rCBF. We estimated κ = 0.51s⁻¹ for the IED-related ANC mechanism around the lesion, which yields $\frac{0.28}{0.51}=0.55$, 14 times smaller than the normal values.

Final Remarks

It is worth noting the foreseeable boost that BOLD modeling will have with the advent of new and optimized sequences in high-field spin-echo fMRI, with the considerably improved ability to measure high-resolution layer-dependent BOLD images and correlates of rCBV and rCBF (47, 70, 84–89). This allows for the construction and estimation of more detailed models of BOLD generation, through understanding of the actual role of arteries, capillaries, and veins in the generation of these observables and the possible biases that the variability of neurovascular/metabolic coupling, CBV, and signal-to-noise ratio (SNR) across layers could introduce. For example, it has been reported that the baseline CBV distribution varies over cortical layers biasing fMRI signal to layers with high CBV values (77). This affects the interpretation of what the contribution of the different vascular compartments to the average low-resolution BOLD response is.

Data Availability Statement

The datasets presented in this article are not readily available because they are from clinical cases and this would jeopardize patient privacy. Requests to access the datasets should be directed to anJpZXJhZGkmI3gwMDA0MDtmaXUuZWR1.

Ethics Statement

The studies involving human participants were reviewed and approved by Western IRB, USA. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

Author Contributions

AS: data analysis and model development. PV-H: data collecting and data analysis. BB: data analysis and clinic case discussion. CD: IED detection and classification. JB-B: data analysis. JR: research design, discussion, and writing. PV-H: data collecting, data analysis, and model development. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Institutes of Health (R56NS094784-01A1).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We appreciate Prof. Pedro A. Valdes-Sosa and Prof. Louis Lemeiux for their useful comments and suggestions during the preparation of this manuscript. Special thanks goes to Professor Jean Gotman and Dr. Nicolas von-Ellenrieder, who substantially helped with the data analysis and interpretation.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fneur.2021.659081/full#supplementary-material

References

1. An D, Fahoum F, Hall J, Olivier A, Gotman J, Dubeau F. Electroencephalography/functional magnetic resonance imaging responses help predict surgical outcome in focal epilepsy. Epilepsia. (2013) 54:2184–94. doi: 10.1111/epi.12434

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Coan AC, Chaudhary UJ, Grouiller F, Campos BM, Perani S, de Ciantis A, et al. EEG-fMRI in the presurgical evaluation of temporal lobe epilepsy. J Neurol Neurosurg Psychiatry. (2016) 87:642–9. doi: 10.1136/jnnp-2015-310401

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Pittau F, Ferri L, Fahoum F, Dubeau F, Gotman J. Contributions of EEG-fMRI to assessing the epileptogenicity of focal cortical dysplasia. Front Comput Neurosci. (2017) 11:8. doi: 10.3389/fncom.2017.00008

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Friston KJ, Holmes AP, Worsley KJ, Poline JB, Frith C, Frackowiak RS. Statistical parametric maps in functional imaging: a general linear approach. Hum Brain Mapping. (1995) 2:189–210. doi: 10.1002/hbm.460020402

CrossRef Full Text | Google Scholar

5. Van Eyndhoven S, Dupont P, Tousseyn S, Vervliet N, Van Paesschen W, Van Huffel S, et al. Augmenting interictal mapping with neurovascular coupling biomarkers by structured factorization of epileptic EEG and fMRI data. NeuroImage. (2021) 228:117652. doi: 10.1016/j.neuroimage.2020.117652

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Lu Y, Grova C, Kobayashi E, Dubeau F, Gotman J. Using voxel-specific hemodynamic response function in EEG-fMRI data analysis: an estimation and detection model. Neuroimage. (2007) 34:195–203. doi: 10.1016/j.neuroimage.2006.08.023

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Grouiller F, Vercueil L, Krainik A, Segebarth C, Kahane P, David O. Characterization of the hemodynamic modes associated with interictal epileptic activity using a deformable model-based analysis of combined EEG and functional MRI recordings. Human Brain Mapp. (2010) 31:1157–73. doi: 10.1002/hbm.20925

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Storti SF, Formaggio E, Bertoldo A, Manganotti P, Fiaschi A, Toffolo GM. Modelling hemodynamic response function in epilepsy. Clin Neurophysiol. (2013) 124:2108–18. doi: 10.1016/j.clinph.2013.05.024

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Proulx S, Safi-Harb M, LeVan P, An D, Watanabe S, Gotman J. Increased sensitivity of fast BOLD fMRI with a subject-specific hemodynamic response function and application to epilepsy. NeuroImage. (2014) 93:59–73. doi: 10.1016/j.neuroimage.2014.02.018

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Bagshaw AP, Aghakhani Y, Bénar C.-G., Kobayashi E, Hawco C, et al. EEG-fMRI of focal epileptic spikes: analysis with multiple haemodynamic functions and comparison with gadolinium-enhanced MR angiograms. Hum Brain Mapp. (2004) 22:179–92. doi: 10.1002/hbm.20024

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Kobayashi E, Bagshaw AP, Grova C, Dubeau F, Gotman J. Negative BOLD responses to epileptic spikes. Hum Brain Mapp. (2006) 27:488–97. doi: 10.1002/hbm.20193

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Rathakrishnan R, Moeller F, Levan P, Dubeau F, Gotman J. BOLD signal changes preceding negative responses in EEG-fMRI in patients with focal epilepsy. Epilepsia. (2010) 51:1837–45. doi: 10.1111/j.1528-1167.2010.02643.x

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Pittau F, Fahoum F, Zelmann R, Dubeau F, Gotman J. Negative BOLD response to interictal epileptic discharges in focal epilepsy. Brain Topogr. (2013) 26:627–40. doi: 10.1007/s10548-013-0302-1

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Salek-Haddadi A, Diehl B, Hamandi K, Merschhemke M, Liston A, Friston K, et al. Hemodynamic correlates of epileptiform discharges: an EEG-fMRI study of 63 patients with focal epilepsy. Brain Res. (2006) 1088:148–66. doi: 10.1016/j.brainres.2006.02.098

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Pesaresi I, Cosottini M, Belmonte G, Maritato P, Mascalchi M, Puglioli M, et al. Reproducibility of BOLD localization of interictal activity in patients with focal epilepsy: intrasession and intersession comparisons. Magma. (2011) 24:285–96. doi: 10.1007/s10334-011-0263-x

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Beers CA, Williams RJ, Gaxiola-Valdez I, Pittman DJ, Kang AT, Aghakhani Y, et al. Patient specific hemodynamic response functions associated with interictal discharges recorded via simultaneous intracranial EEG-fMRI. Hum Brain Mapp. (2015) 36:5252–64. doi: 10.1002/hbm.23008

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Riera JJ, Watanabe J, Kasuki I, Naoki M, Aubert E, Ozaki T, et al. A state-space model of the hemodynamic approach: nonlinear filtering of BOLD signals. NeuroImage. (2004) 21:547–67. doi: 10.1016/j.neuroimage.2003.09.052

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Makni S, Ciuciu P, Idier J, Poline JB. Joint detection-estimation of brain activity in functional MRI: a multichannel deconvolution solution. IEEE Trans Signal Process. (2005) 53:3488–502. doi: 10.1109/TSP.2005.853303

CrossRef Full Text | Google Scholar

19. Makni A, Idier J, Vincent T, Thirion B, Dehaene-Lambertz G, Ciuciu P. A fully Bayesian approach to the parcel-based detection-estimation of brain activity in fMRI. NeuroImage. (2008) 41:941–69. doi: 10.1016/j.neuroimage.2008.02.017

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Baraldi P, Manginelli AA, Maieron M, Liberati D, Porro CA. An ARX model-based approach to trial by trial identification of fMRI-BOLD responses. Neuroimage. (2007) 37:189–201. doi: 10.1016/j.neuroimage.2007.02.045

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Vincent T, Risser L, Ciuciu P. Spatially adaptive mixture modeling for analysis of fMRI time series. IEEE Trans Medical Imaging. (2010) 29:1059–74. doi: 10.1109/TMI.2010.2042064

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Chaari L, Vincent T, Forbes F, Dojat M, Ciuciu P. Fast joint detection-estimation of evoked brain activity in event-related fMRI using a variational approach. IEEE Trans Medical Imaging. (2013) 32:821–37. doi: 10.1109/TMI.2012.2225636

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Cherkaoui H, Moreau T, Halimi A, Leroy C, Ciuciu P. Multivariate semi-blind deconvolution of fMRI time series. Neuroimage. (2021) 241:118418. doi: 10.1016/j.neuroimage.2021.118418

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Lemieux L, Salek-Haddadi A, Josephs O, Allen P, Toms N, Scott C, et al. Event-related fMRI with simultaneous and continuous EEG: description of the method and initial case report. Neuroimage. (2001) 14:780–7. doi: 10.1006/nimg.2001.0853

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Neckelmann D, Amzica F, Steriade M. Changes in neuronal conductance during different components of cortically generated spike-wave seizures. Neuroscience. (2000) 96:475–85. doi: 10.1016/S0306-4522(99)00571-0

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Pollen DA. Intracellular studies of cortical neurons during thalamic induced wave and spike. Electroencephalogr Clin Neurophysiol. (1964) 17:398–404. doi: 10.1016/0013-4694(64)90163-4

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Maggioni E, Zucca C, Reni G, Cerutti S, Triulzi FM, Bianchi AM, et al. Investigation of the electrophysiological correlates of negative BOLD response during intermittent photic stimulation: an EEG-fMRI study. Hum Brain Mapp. (2016) 37:2247–62. doi: 10.1002/hbm.23170

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Pasley BN, Inglis BA, Freeman RD. Analysis of oxygen metabolism implies a neural origin for the negative BOLD response in human visual cortex. Neuroimage. (2007) 36:269–76. doi: 10.1016/j.neuroimage.2006.09.015

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Schäfer K, Blankenburg F, Kupers R, Grüner JM, Law I, Lauritzen M, et al. Negative BOLD signal changes in ipsilateral primary somatosensory cortex are associated with perfusion decreases and behavioral evidence for functional inhibition. Neuroimage. (2012) 59:3119–27. doi: 10.1016/j.neuroimage.2011.11.085

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Shmuel A, Yacoub E, Pfeuffer J, Van de Moortele PF, Adriany G, Hu X, et al. Sustained negative BOLD, blood flow and oxygen consumption response and its coupling to the positive response in the human brain. Neuron. (2002) 36:1195–210. doi: 10.1016/S0896-6273(02)01061-9

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Smith AT, Williams AL, Singh KD. Negative BOLD in the visual cortex: evidence against blood stealing. Hum Brain Mapp. (2004) 21:213–20. doi: 10.1002/hbm.20017

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Wade AR. The negative BOLD signal unmasked. Neuron. (2002) 36:993–5. doi: 10.1016/S0896-6273(02)01138-8

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Van Paesschen W. Ictal SPECT. Epilepsia. (2004) 45:35–40. doi: 10.1111/j.0013-9580.2004.04008.x

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Sarikaya I. PET studies in epilepsy. Am J Nucl Med Mol Imaging. (2015) 5:416–30.

Google Scholar

35. Song Y, Torres R, Garcia S, Frometa Y, Bae J, Deshmukh A, et al. Dysfunction of neuro-vascular/metabolic coupling in chronic focal epilepsy. IEEE Trans Biomed Eng. (2016) 63:97–110. doi: 10.1109/TBME.2015.2461496

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Fahoum F, Zelmann R, Tyvaert L, Dubeau F, Gotman J. Epileptic discharges affect the default mode network - fMRI and intracerebral EEG evidence. PLoS One. (2013) 8:e68038. doi: 10.1371/journal.pone.0068038

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Laufs H, Hamandi K, Salek-Haddadi A, Kleinschmidt AK, Duncan JS, Lemieux L. Temporal lobe interictal epileptic discharges affect cerebral activity in “default mode” brain regions. Hum Brain Mapp. (2007) 28:1023–32. doi: 10.1002/hbm.20323

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Harel N, Lee SP, Nagaoka T, Kim DS, Kim SG. Origin of negative blood oxygenation level-dependent fMRI signals. J Cereb Blood Flow Metab. (2002) 22:908–17. doi: 10.1097/00004647-200208000-00002

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Suarez A, Valdes-Hernandez PA, Moshkforoush A, Tsoukias N, Riera J. Arterial blood stealing as a mechanism of negative BOLD response: from the Steady-flow with nonlinear phase separation to a Windkessel-based model. J Theor Biol. (2021) 529. doi: 10.1016/j.jtbi.2021.110856

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Hu D, Huang L. Negative hemodynamic response in the cortex: evidence opposing neuronal deactivation revealed via optical imaging and electrophysiological recording. J Neurophysiol. (2015) 114:2152–61. doi: 10.1152/jn.00246.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Ma Z, Cao P, Sun P, Lu Z, Li L, Chen Y, et al. Negative hemodynamic response without neuronal inhibition investigated by combining optical imaging and electrophysiological recording. Neurosci Lett. (2017) 637:161–7. doi: 10.1016/j.neulet.2016.11.029

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Mandeville JB, Marota JJ, Ayata C, Zaharchuk G, Moskowitz MA, Rosen BR, et al. Evidence of a cerebrovascular postarteriole Windkessel with delayed compliance. J Cereb Blood Flow Metab. (1999) 19:679–89. doi: 10.1097/00004647-199906000-00012

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Buxton RB, Wong EC, Frank LR. Dynamics of blood flow and oxygenation changes during brain activation: the Balloon model. MRM. (1998) 39:855–64. doi: 10.1002/mrm.1910390602

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Havlicek M, Ivanov D, Roebroeck A, Uludag K. Determining excitatory and inhibitory neuronal activity from multimodal fMRI data using a generative hemodynamic model. Front Neurosci. (2017) 11:616. doi: 10.3389/fnins.2017.00616

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Riera J, Bosch J, Yamashita O, Kawashima R, Sadato N, Okada T, et al. fMRI activation maps based on the NN-ARx model. NeuroImage. (2004) 23:680–97. doi: 10.1016/j.neuroimage.2004.06.039

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Friston K. Statistical parametric mapping. In: Statistical Parametric Mapping: The Analysis of Functional Brain Images. London: Academia Press; Elsevier (2007). p. 10–31.

Google Scholar

47. Goense J, Merkle H, Logothetis NK. High-resolution fMRI reveals laminar differences in neurovascular coupling between positive and negative BOLD responses. Neuron. (2012) 76:629–39. doi: 10.1016/j.neuron.2012.09.019

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Ashburner J, Friston KJ. Unified segmentation. Neuroimage. (2005) 26:839–51. doi: 10.1016/j.neuroimage.2005.02.018

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Allen PJ, Josephs O, Turner R. A method for removing imaging artifact from continuous EEG recorded during functional MRI. Neuroimage. (2000) 12:230–9. doi: 10.1006/nimg.2000.0599

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Allen PJ, Polizzi G, Krakow K, Fish DR, Lemieux L. Identification of EEG events in the MR scanner: the problem of pulse artifact and a method for its subtraction. Neuroimage. (1998) 8:229–39. doi: 10.1006/nimg.1998.0361

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Bell AJ, Sejnowski TJ. An information-maximization approach to blind separation and blind deconvolution. Neural Comput. (1995) 7:1129–59. doi: 10.1162/neco.1995.7.6.1129

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Makeig S, Jung TP, Bell AJ, Ghahremani D, Sejnowski TJ. Blind separation of auditory event-related brain responses into independent components. Proc Natl Acad Sci USA. (1997) 94:10979–84. doi: 10.1073/pnas.94.20.10979

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Havlicek M, Roebroeck A, Friston K, Gardumi A, Ivanov D, Uludag K. Physiologically informed dynamic causal modeling of fMRI data. Neuroimage. (2015) 122:355–72. doi: 10.1016/j.neuroimage.2015.07.078

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Zheng Y, Mayhew J. A time-invariant visco-elastic Windkessel model relating blood flow and blood volume. Neuroimage. (2009) 47:1371–80. doi: 10.1016/j.neuroimage.2009.04.022

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Friston KJJ, Mechelli A, Turner R, Price CJJ. Nonlinear responses in fMRI: the Balloon model, Volterra kernels, and other hemodynamics. Neuroimage. (2000) 12:466–77. doi: 10.1006/nimg.2000.0630

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Zheng Y, Martindale J, Johnston D, Jones M, Berwick J, Mayhew J. A model of the hemodynamic response and oxygen delivery to brain. Neuroimage. (2002) 16:617–37. doi: 10.1006/nimg.2002.1078

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Glover GH. Deconvolution of impulse response in event-related BOLD fMRI. Neuroimage. (1999) 9:416–29. doi: 10.1006/nimg.1998.0419

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Chen C-C, Tyler CW. Spectral analysis of fMRI signal noise. In: Onozuka M, Yen C-T, editors. Novel Trends in Brain Science: Brain Imaging, Learning and Memory, Stress and Fear, and Pain. Tokyo: Springer Japan (2008). p. 63–76.

Google Scholar

59. Friston KJJ. Bayesian estimation of dynamical systems: an application to fMRI. Neuroimage. (2002) 16:513–30. doi: 10.1006/nimg.2001.1044

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Furey TS, Cristianini N, Duffy N, Bednarski DW, Schummer M, Haussler D. Support vector machine classification and validation of cancer tissue samples using microarray expression data. Bioinformatics. (2000) 16:906–14. doi: 10.1093/bioinformatics/16.10.906

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Mathur A, Foody GM. Multiclass and binary SVM classification: implications for training and classification users. IEEE Geosci Remote Sensing Lett. (2008) 5:241–5. doi: 10.1109/LGRS.2008.915597

CrossRef Full Text | Google Scholar

62. Cervantes J, Garcia-Lamont F, Rodríguez-Mazahua L, Lopez A. A comprehensive survey on support vector machine classification: applications, challenges and trends. Neurocomputing. (2020) 408:189–215. doi: 10.1016/j.neucom.2019.10.118

CrossRef Full Text | Google Scholar

63. Ayala GF, Dichter M, Gumnit RJ, Matsumoto H, Spencer WA. Genesis of epileptic interictal spikes. New knowledge of cortical feedback systems suggests a neurophysiological explanation of brief paroxysms. Brain Res. (1973) 52:1–17. doi: 10.1016/0006-8993(73)90647-1

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Deco G, Jirsa VK, McIntosh AR. Emerging concepts for the dynamical organization of resting-state activity in the brain. Nat Rev Neurosci. (2011) 12:43–56. doi: 10.1038/nrn2961

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Deco G, Jirsa VK. Ongoing cortical activity at rest: criticality, multistability, and ghost attractors. J Neurosci. (2012) 32:3366–75. doi: 10.1523/JNEUROSCI.2523-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Ghosh A, Rho Y, McIntosh AR, Kötter R, Jirsa VK. Noise during rest enables the exploration of the brain's dynamic repertoire. PLoS Comput Biol. (2008) 4:e1000196. doi: 10.1371/journal.pcbi.1000196

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Haken H. Synergetics: An Introduction: Nonequilibrium Phase Transitions and Self-Organization in Physics, Chemistry, and Biology, 3rd ed, Springer Series in Synergetics. New York, NY: Springer-Verlag (1983).

Google Scholar

68. Bittar RG, Andermann F, Olivier A, Dubeau F, Dumoulin SO, Pike GB, et al. Interictal spikes increase cerebral glucose metabolism and blood flow: a PET study. Epilepsia. (1999) 40:170–8. doi: 10.1111/j.1528-1157.1999.tb02071.x

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Boas DA, Jones SR, Devor A, Huppert TJ, Dale AM. A vascular anatomical network model of the spatio-temporal response to brain activation. Neuroimage. (2008) 40:1116–29. doi: 10.1016/j.neuroimage.2007.12.061

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Bandettini PA. The BOLD plot thickens: sign- and layer-dependent hemodynamic changes with activation. Neuron. (2012) 76:468–9. doi: 10.1016/j.neuron.2012.10.026

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Buxton RB, Uludag K, Dubowitz DJ, Liu TT. Modeling the hemodynamic response to brain activation. Neuroimage. (2004) 23:S220–33. doi: 10.1016/j.neuroimage.2004.07.013

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Shmuel A, Augath M, Oeltermann A, Logothetis NK. Negative functional MRI response correlates with decreases in neuronal activity in monkey visual area V1. Nat Neurosci. (2006) 9:569–77. doi: 10.1038/nn1675

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Geyer S, Schormann T, Mohlberg H, Zilles K. Areas 3a, 3b, and 1 of human primary somatosensory cortex. Neuroimage. (2000) 11:684–96. doi: 10.1006/nimg.2000.0548

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Cappe C, Rouiller E, Barone P. Cortical and thalamic pathways for multisensory and sensorimotor interplay. In: The Neural Bases of Multisensory Processes. Boca Raton, FL: CRC Press/Taylor & Francis (2011). p. 15–30.

Google Scholar

75. Hagmann P, Kurant M, Gigandet X, Thiran P, Wedeen VJ, Meuli R, et al. Mapping human whole-brain structural networks with diffusion MRI. PLoS ONE. (2007) 2:e597. doi: 10.1371/journal.pone.0000597

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Joshi AA, Joshi SH, Dinov I, Shattuck DW, Leahy RM, Toga AW. Anatomical structural network analysis of human brain using partial correlations of gray matter volumes. In: 2010 IEEE International Symposium on Biomedical Imaging: From Nano to Macro. Rotterdam: IEEE (2010). p. 844–7.

Google Scholar

77. Uludag K, Blinder P. Linking brain vascular physiology to hemodynamic response in ultra-high field MRI. Neuroimage. (2018) 168:279–95. doi: 10.1016/j.neuroimage.2017.02.063

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Fahoum F, Lopes R, Pittau F, Dubeau F, Gotman J. Widespread epileptic networks in focal epilepsies: EEG-fMRI study. Epilepsia. (2012) 53:1618–27. doi: 10.1111/j.1528-1167.2012.03533.x

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Norden AD, Blumenfeld H. The role of subcortical structures in human epilepsy. Epilepsy Behav. (2002) 3:219–31. doi: 10.1016/S1525-5050(02)00029-X

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Struck AF, Boly M, Hwang G, Nair V, Mathis J, Nencka A, et al. Regional and global resting-state functional MR connectivity in temporal lobe epilepsy: results from the Epilepsy Connectome Project. Epilepsy Behav. (2021) 117:107841. doi: 10.1016/j.yebeh.2021.107841

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Schridde U, Khubchandani M, Motelow JE, Sanganahalli BG, Hyder F, Blumenfeld H. Negative BOLD with large increases in neuronal activity. Cereb Cortex. (2008) 18:1814–27. doi: 10.1093/cercor/bhm208

PubMed Abstract | CrossRef Full Text | Google Scholar

82. Stefanovic B, Warnking JM, Kobayashi E, Bagshaw AP, Hawco C, Dubeau F, et al. Hemodynamic and metabolic responses to activation, deactivation and epileptic discharges. Neuroimage. (2005) 28:205–15. doi: 10.1016/j.neuroimage.2005.05.038

PubMed Abstract | CrossRef Full Text | Google Scholar

83. Gallagher A, Madan N, Stemmer-Rachamimov A, Thiele EA. Progressive calcified tuber in a young male with tuberous sclerosis complex. Dev Med Child Neurol. (2010) 52:1062–65. doi: 10.1111/j.1469-8749.2010.03792.x

PubMed Abstract | CrossRef Full Text | Google Scholar

84. Feinberg DA, Yacoub E. The rapid development of high speed, resolution and precision in fMRI. Neuroimage. (2012) 62:720–5. doi: 10.1016/j.neuroimage.2012.01.049

PubMed Abstract | CrossRef Full Text | Google Scholar

85. Harel N, Bolan PJ, Turner R, Ugurbil K, Yacoub E. Recent advances in high-resolution MR application and its implications for neurovascular coupling research. Front Neuroenergetics. (2010) 2:130. doi: 10.3389/fnene.2010.00130

PubMed Abstract | CrossRef Full Text | Google Scholar

86. Huber L, Goense J, Kennerley AJ, Ivanov D, Krieger SN, Lepsien J, et al. Investigation of the neurovascular coupling in positive and negative BOLD responses in human brain at 7T. Neuroimage. (2014) 97:349–62. doi: 10.1016/j.neuroimage.2014.04.022

PubMed Abstract | CrossRef Full Text | Google Scholar

87. Huber L, Goense J, Kennerley AJ, Trampel R, Guidi M, Reimer E, et al. Cortical lamina-dependent blood volume changes in human brain at 7T. Neuroimage. (2015) 107:23–33. doi: 10.1016/j.neuroimage.2014.11.046

PubMed Abstract | CrossRef Full Text | Google Scholar

88. Polimeni JR, Uludag K. Neuroimaging with ultra-high field MRI: present and future. Neuroimage. (2018) 168:1–6. doi: 10.1016/j.neuroimage.2018.01.072

PubMed Abstract | CrossRef Full Text | Google Scholar

89. Goense J, Bohraus Y, Logothetis NK. fMRI at High Spatial Resolution: Implications for BOLD-Models. Front Comput Neurosci. (2016) 10:66. doi: 10.3389/fncom.2016.00066

PubMed Abstract | CrossRef Full Text | Google Scholar