Noninvasive Epileptogenic Signal Direction Determination by Effective Connectivity of Resting State Functional MRI - Independent of EEG

1. Introduction¹

The most effective and only known curative treatment for drug resistant epilepsy (DRE) is surgery to remove the seizure focus and or interrupt the epileptogenic network (Luders et al., 2006). The primary determinant of surgical candidacy and success is accurate localization of the seizure onset zone (SOZ) (West et al., 2019). However, current noninvasive SOZ localization methods still often depend on confirmation by intracranial electroencephalograph (iEEG), which is expensive, carry risks, and still only leads to 40-80% seizure-freedom when a seizure focus is “identified” (Tonini et al., 2004).

Improvement in noninvasive SOZ source localization over more standard methods has been demonstrated by static connectivity from resting state fMRI (rs-fMRI) via independent component analysis (ICA) (Chakraborty et al., 2020). Notably, rs-fMRI ICA-derived SOZ findings have not only been associated with iEEG, but also with increased surgical candidacy. Additionally, improved Engel outcomes, have been associated with resolution of post-operative rs-fMRI SOZ ICA networks (Boerwinkle et al., 2019b; Boerwinkle et al., 2020; Boerwinkle et al., 2017). However, after evaluating over 2000 rs-fMRI of individual studies with DRE using ICA for SOZ localization, a major weakness remains: this static network measure identifies more than one plausible SOZ candidate in at least 40% of subjects (VB). Thus, there remains dependency on stereotactic iEEG to distinguish the best surgical target a concordance of data from rs-fMRI and other noninvasive methods.

The fundamental limitation of static rs-fMRI ICA is overcome to some extent by employing simultaneous EEG-fMRI, though this method only identifies the SOZ in a small proportion of cases, alone (Vaudano et al., 2021). The interictal epileptogenic discharges (IED) detected by EEG inform the time point at which pre-to-post blood oxygenation level dependent (BOLD) imaging maximal regions of interest (IED-BOLD ROI) are identified. A possible advance on this method, subsequent dynamic causal modeling (DCM) of the IED-BOLD ROI, has differentiated the epileptogenic driver from propagation regions in four cases as determined by good surgical outcomes (Engel class I and II) (Engel, 1993; Vaudano et al., 2021). However, this limited series highlighted the dependency of EEG-fMRI on requires additional specialized equipment, EEG technician staff, and epileptiform events to occur during the scan; all of which are cumbersome features limiting clinical utility.

To overcome these limitations, we designed a computational approach, Directionality, which employs cross-spectral DCM from rs-fMRI alone, without electrophysiological information, known as resting state effective connectivity (RSEC) helps to differentiate the driving SOZ region from propagation regions. These areas can be initially identified from the static rs-fMRI analysis. Directionality assumes SOZ are: First, (1) generators of excitatory (positive) signal toward regions of propagation; and secondly (2) receivers of inhibitory (negative) activity from the propagation node(s).

To determine if Directionality accurately distinguishes between SOZ and propagation zones (pZ), we selected a DRE population with homogenous and established SOZ localization, as well as prior surgically validated static rs-fMRI SOZ and pZ determined by the static rs-fMRI connectivity measure, SearchLight (SL) (Boerwinkle et al., 2018a). DRE from our patients with hypothalamic hamartoma (HH-DRE) fills these criteria because (1) it is well-established that the HH is the primary seizure driver in HH-DRE; (2) the HH epileptogenic network dynamic follow the Papez circuit, as verified by prior iEEG, rs-fMRI by ICA and partial correlation (Boerwinkle et al., 2016) and SL (Boerwinkle et al., 2018a), and EEG-fMRI by DCM (Murta et al., 2012; Usami et al., 2016). Notably, the same pathway of seizure progression from the HH to the rest of the brain from iEEG was found by rs-fMRI static measures alone in HH, but only on a group level (Boerwinkle et al., 2016) inferring the potential success of SOZ-identifying information from rs-fMRI Directionality on an individual basis.

In this study, we tested our hypotheses of signal direction by Directionality to distinguish the known HH’s SOZ from an ROI of propagation (the pZ) in children with HH-DRE who underwent pre-operative rs-fMRI SL-guided laser interstitial thermal ablation therapy (LITT) and had Engel I surgical outcomes. Further, we hypothesized that higher strength of signal from the identified region with negative signal (the pZ) and higher excitation from the region with positive signal (the SOZ) will correlate with seizure outcome. Also, because length of time of DRE is associated with decreased improvement with surgery (Tonini et al., 2004), we hypothesized that age will be associated with signal strength between these regions.

2. Methods

2.2 MRI Data acquisition and processing

2.2.1 MRI acquisition

Images were acquired on a 3 Tesla MRI (Ingenuity, Philips Medical Systems, Best, Netherlands) equipped with a 32-channel head coil. Rs-fMRI parameters included TR (repetition time) 2000 ms, TE (echo time) 30 ms, matrix size 80×80, flip angle 80, number of slices 46, slice thickness 3.4 mm with no gap, in plane resolution 3×3 mm, inter-leaved acquisition, and number of total volumes 600, in two 10 minute runs, totaling 20 minutes. For anatomical reference, a T1-weighted turbo field echo whole-brain sequence was obtained with the following parameters: TR 9 ms, TE 4 ms, flip angle 8, slice thickness 0.9 mm, and in-plane resolution 0.9×0.9 mm.

2.2.3 MRI preprocessing

Using the program Statistical Parametric Mapping version 12 (Friston et al., 1994) (SPM12; www.fil.ion.ucl.ac.uk/spm, Wellcome Trust Centre for Neuroimaging, London, UK), in Matlab 2019b, the T1 image was resampled to 1×1×1 mm voxels, the origin was set at the anterior commissure and, if needed, minor adjustments were made to correct orientation (head-tilt). The T1 image was segmented using the Computational Anatomy Toolbox version 12.7 (CAT12) (Gaser and Dahnke, 2016) and the Automatic Anatomical Labelling Atlas 3 (AAL3) (Rolls et al., 2020) was registered into subject T1 space and masks for each of the hippocampi were generated.

For the fMRI data, the same pre-processing steps as prior work (Boerwinkle et al., 2019b; Boerwinkle et al., 2017) were applied by high-pass filtering the data to remove ultra-low-frequency non-neural artifacts and extract the grey matter voxel time-course while removing voxels in the cerebrospinal fluid, and correcting for subject movement. Following ICA-based denoising (manual expert classification with removal of noise-based components in epilepsy (Boerwinkle et al., 2019b; Boerwinkle et al., 2017), CSF, and motion regression), the denoised functional data were re-aligned in SPM to register both of each subjects’ 10-minute rs-fMRI runs to each other. The mean functional image was segmented with unified segmentation. Using the bias-corrected mean functional image, the origin was set and orientation adjustments were made as described above for the T1. Functional images were co-registered (estimate) to T1 space in SPM using the bias-corrected mean functional image, with visual inspection.

2.2.4 ROI Selection

Network ROI locations selected were: (1) the SOZ located by prior clinical SL and visibly ablated as shown in the immediate post-operative diffusion weighted sequencing, and (2) the region with the largest intracerebral cluster of HH-connectivity by SL, selecting the highest intra-cluster connectivity by SL (patient example in Figure 1, all patients in Supplementary Table 1). The choice and size of the SL-guided SOZs were determined by the constraints of (1) the HH size, (2) the region within the ablated portion of the HH that preoperatively had the largest and highest activation, and (3) highest correlation to the regions outside the HH in pZ pattern demonstrated by SL. The size of the pZ was primarily selected for regional grey matter thickness within this peak-appearing connectivity region, or otherwise, region of cohesive and high connectivity determined by expert review and known to be a frequent region of propagation by prior studies, and primarily the individual’s SL results (Boerwinkle et al., 2018a; Boerwinkle et al., 2016; Murta et al., 2012; Usami et al., 2016). The pZ was either drawn manually or an AAL3 hippocampus mask. Time series extraction from the SOZ and pZ were performed through a general linear model (GLM) of the rs-fMRI data. Given the prior noise correction, no additional noise correction regressors were used. Both rs-fMRI runs were entered as the same session. SPM concatenate (Casanova et al., 2007) was used to specify the boundary between the runs to correct for intensity differences between sessions. From the resulting SPM, a time series for each SOZ and pZ were extracted using the first eigenvector, as prior (Sussman et al., 2021).

• Download figure
• Open in new tab

Figure 1.

Example of ROI of the SOZ and pZ selection. Each set of six images is divided into two rows of sagittal, coronal, and axial images. A. Pre-operative T1W wherein row 1 is a left hippocampal pZ (red circle), and in row 2 is the SOZ within the HH, which is one voxel large. B-D. Post-operative post-contrast T1W, pre-operative T1W, and post-operative diffusion weighed images, respectively showing in row 1 same SOZ as in A, and row 2 includes the SL connectivity to the SOZ. These views allow visualization of the regions of the SOZ destroyed as either the post-contrast or the diffusion were positive in SOZ location and with connectivity to the pZ ROI selected. Abbreviations: ROI: region of interest; pZ: propagation zone; SOZ: seizure onset zone; SL: SearchLight.

2.2.5 DCM

2.2.5.1 Model Estimation

In overview, cross-spectral DCM (Friston et al., 2014a) was used to model the regional time course signals, in a fully-connected model, similar to previously described (Sussman et al., 2021). Cross-spectral DCM estimates parameters of auto- and cross-spectrum through multivariate auto-regression models of BOLD data generating estimated spectrums. As such, a two-node model was specified per case and steady-state spectral amplitude and phase representations of each node’s (ROI) first principle component activity were obtained through Fourier transform. Variational Bayesian inversion was used to fit the differential equation connectivity model (Friston et al., 2003; Friston et al., 2014b). Inversion involves fitting the DCM to maximize the likelihood of the model under prior parameter specifications. Model-fitting estimates parameters that describe the amplitude by frequency-spectral representation for each region. Each region’s local spectrum was modeled as a power law distribution with an amplitude and scale, the latter indicating the frequency by amplitude slope. Effective connectivity was derived through estimating the same parameters through frequency cross-spectrum between regions.

2.2.5.2 Single Subject Model Comparison

Our aim of identifying connections in this HH network was achieved by a Bayesian model reduction (BMR) and averaging approach (BMA). Specifically, we estimated the probability that every single model generated the data through its Bayesian furnished model evidence of the fully-connected model (BMR). The model evidence was then used to weight the probability of each connectivity parameter in each fit model. These parameters were then averaged across all models to determine the statistically optimal estimate of local neural and regional connectivity effects (BMA). This approach allowed us to furnish a statistical reliability measure of estimated connectivity parameters for each tested model within a single subject, circumventing the need to compare models across subjects. A threshold of 0.9 was used for the posterior probability of each parameter for further analysis.

2.2.5.3 Group Model Estimation and Comparison

In addition to single-subject model comparison, we planned a group Parametric Empirical Bayes (PEB) analysis to perform group-based model reduction and comparison to identify mean group effects as well as investigate the group effects of Engel class outcome and age on parameter estimate size. For this, the fully-connected individual DCMs were estimated, but then exhaustive Bayesian model reduction, comparison, and averaging were performed through a PEB paradigm (Friston et al., 2016; Friston et al., 2015; Zeidman et al., 2019). Exhaustive BMR was used because it is a data-driven strategy in an otherwise hypothesis driven technique. In the model, Engel class 1a outcome were coded as 1 and Engel class 1b outcome were coded as -1. The covariates were mean-centered and entered in the following order: Engel class outcome, age. Covariate effects of interest (posterior probability > 0.95) were entered into a leave-one-out (LOO) cross-validation analysis to examine their predictive ability. For discussion purposes, we focused on parameters with a posterior probability greater than 0.95 (free energy).

3. Results

3.2 Comparisons with individually estimated and reduced models

Results from the individual DCM estimation, model reduction, and Bayesian model averaging for an example subject are shown in Figure 2. Individual model results for all subjects can be found in Supplementary Table 1. When considering the SOZ effective connectivity, 3.6% had no significant connection (Table 2). Of the remainder with significant connections, 96% and 3.7% were in the expected and reversed directions, respectively (Table 2). Positive and negative polarity signals were equivalent as expected between the SOZ and pZ (Table 2). By connection, considering both positive and negative, Directionality had sensitivity of 96%, accuracy of 93%, and positive predictive value of 96%. The observed agreement between Directionality and surgical outcome was 92.9% (N=52/56; 95% CI: 0.83-0.98). The PABAK estimate was 0.86 (95% CI: 0.65-0.96), indicating substantial agreement between Directionality and surgical outcome.

• Download figure
• Open in new tab

Figure 2.

Individual DCM outputs. These outputs were calculated for each patient. A. Parameter estimates and their posterior probabilities of the fully specified DCM A matrix (intrinsic connectivity) after inversion using cross-spectral DCM. Target region is listed on the x-axis; connections from the SOZ (HH) are in blue and connections from the pZ are in grey. The left graph is the effect sizes, self-connections are shown in the log scaling value used to stabilize the model (conversion to Hz is -0.5 * exp(A)). The right bar chart shows the posterior probabilities of each estimated parameter. B. Adjacency matrix of the A matrix parameter estimates after optimization (exhaustive BMR followed by BMA). This matrix is thresholded to only include connections with a posterior probability greater than 0.9. Self-connections have been converted from the log-scaling value to Hz. C. The values from panel C shown as a diagram superimposed on the patient’s coronal T1W image. Abbreviations: DCM: dynamic causal modeling; ROI: region of interest; pZ: propagation zone; SOZ: seizure onset zone; SL: SearchLight; HH: hypothalamic hamartoma.

View this table:

• View inline
• View popup
• Download powerpoint

Table 2.

Directionality Results and Signal Characterization

3.3 Parametric Empirical Bayes Group Analysis Results

The BMA group model is shown in Figure 3A. The connection from the SOZ to the pZ was excitatory and the connection from the pZ→SOZ was inhibitory. The SOZ self-connection was also present (inhibitory). The pZ self-connection was pruned during BMR.

• Download figure
• Open in new tab

Figure 3.

A. Group baseline connections between the SOZ and pZ thresholded at 0.95 posterior probability after exhaustive Bayesian Model Reduction (BMR) and Bayesian Model Averaging (BMA). The location of the pZ varied by subject. Red arrows and values are excitatory modulation, blue arrows and values are inhibitory modulation. B. Group connections associated with Engel outcome and age thresholded at 0.95 posterior probability. Blue arrow/value indicates a negative correlation and a red arrow/value indicates a positive correlation. B1. Engel class Ia outcomes were coded with a higher value than Engel class Ib outcomes, thus, the blue arrow in Panel A indicates that patients with better one year outcomes (Engel class Ia) showed more inhibitory modulation from the pZ to SOZ. B2. Older patients show less inhibition from the pZ to SOZ and less excitation from the SOZ to pZ, thus, overall, absolute strength of these parameter effect sizes decreased with age

Regarding behavior correlations, Engel class outcome was associated with signal strength from the pZ to the SOZ (Figure 3B). Patients with a better Engel class outcome (1a) showed greater inhibitory signal from the ROI to the SOZ (Figure 3B1). However, this difference, however, was not significant with LOO cross-validation (r(26)=-0.07, P=0.65; Supplementary Figure 1A). Age was also associated with between-node connections; age was positively correlated with signal from the pZ→SOZ and negatively correlated with signal from the SOZ→pZ (Figure 3B2). Neither of these correlations were significant with LOO cross-validation r(26)=0.08, P=0.33 and r(26)=0.21, P=0.14, respectively (Supplementary Figure 1B-C). Self-connections were not related to Engel outcome or age.

4. Discussion

This is the first reported DCM analysis of rs-fMRI, performed independent of electrophysiological data, to demonstrate the capacity of RSEC alone to differentiate SOZ from the region of propagation with high accuracy, sensitivity, and positive predictive value. This is also the first study to demonstrate that cross-spectral DCM of rs-fMRI yields the excitatory signal direction from SOZ to a region of propagation and inhibitory signal from the region of propagation back to the SOZ, on both in the individual and group level. Before application to the more heterogeneously localized SOZ of the general DRE population, it was necessary to perform this analysis in such a homogenous and well established SOZ location of the HH to understand the potential of this tool for narrowing ROI for SOZ location determination. Importantly, polarity of SOZ and pZ connections were consistent in both individual and group analyses, increasing confidence in individual application, similar to prior RSFC (Boerwinkle et al., 2019a; Boerwinkle et al., 2018a; Boerwinkle et al., 2017; Boerwinkle et al., 2016) and RSEC work (Vaudano et al., 2021).

Our hypothesis that, from rs-fMRI alone, intrinsic modulation from the SOZ is excitatory and modulation from the pZ is inhibitory is supported by our results. Unlike previous studies, we did not time lock or model around ictal onset or known periods of interictal discharge (Daunizeau et al., 2012; Hamandi et al., 2008; Klamer et al., 2018; Klamer et al., 2015; Murta et al., 2012; Vaudano et al., 2013; Vaudano et al., 2009; Vaudano et al., 2021; Warren et al., 2019). Instead, we assumed, that in DRE, markers of epileptic signal are present in resting-state networks regardless of known or unknown ictal and interictal activity at time of data capture and without time-locking to ictal/interictal activity. This has been supported in previous studies using RSFC to reliably identify the (static) spatial locations of epileptogenic networks (but not the direction of signal propagation) between regions in these networks without time-locking to epileptogenic activity (Boerwinkle et al., 2018a; Boerwinkle et al., 2017; Boerwinkle et al., 2018b; Boerwinkle et al., 2016). We expanded on previous findings to show that information about direction of epileptogenic signal is also present in event-agnostic rs-fMRI signal.

While it is possible to use methodologies to know whether patients were seizing in the scanner or would have shown concurrent epileptiform activity, the consistency of the results herein is a strength that implies that, in this population, such knowledge may not be necessary in order to model seizure propagation direction but requires further study in a larger population of these patients who had poor outcomes and see if the model remains consistent. In fact, necessitating EEG-fMRI and capture of epileptiform event in scanner can reduce the overall clinical yield of fMRI-based EC techniques in epilepsy network mapping. Using EEG-fMRI to guide DCM, Vaudano et al. (2021) found that all 4 patients whose DCM indicated SOZ was concordant with the clinical SOZ and surgically-targeted had good clinical outcomes. However, this subsample represented only a 14.3% (4/28) surgically-validated clinical yield of all patients who underwent EEG-fMRI and would likely benefit from an EC analysis, compared to a surgically-validated clinical yield of 81.3% (26/32) in the current study (see Supplementary Table 6), a Fischer exact test found this difference statistically significant (P<0.001). A broad summary of methodological differences between the current study and prior studies is detailed in Table 3.

Our current results show some differences in terms of polarity from prior studies. For example, we show that the connection from pZ→SOZ is overwhelmingly negative; however, pZ→SOZ connections in prior event-based DCM studies in have been mostly reported positive (Hamandi et al., 2008; Murta et al., 2012; Vaudano et al., 2013) with one study showing a mix (Vaudano et al., 2013). Thus, our results show greater consistency of A matrix connection polarities. This difference may be due to the number of nodes modeled, differences in ROI specification, patient effects (all of our patients had Engel class I outcomes), or potentially more likely – our model did not include direct inputs or modulatory influences during captured explicit epileptiform activity and used cross-spectra (second-order statistic) instead of time-varying fluctuations (time series). As such, we also assume that the state of the relationship between nodes in the current study is modelling the basal state of the relationship between SOZ and pZ without placing temporal emphasis on ictal or inter-ictal events. This may be different from their relationship during these events, thus contributing to differences seen in our analysis. However, the differences could also be related to volatility of nodal relationship during epileptiform activity, as nodal strength and polarity inconsistencies in prior EEG-fMRI are prevalent (David et al., 2008; Murta et al., 2012; Vaudano et al., 2013), whereas Directionality was highly consistent.

Also consistent for Directionality was network architecture. Bi-directional connections between SOZ and pZ had overwhelmingly high probability, which is not dissimilar previous time-locked DCM studies that compared models with and without bidirectional connections along propagation pathways (Murta et al., 2012; Warren et al., 2019). Thus, our results imply that both SOZ→pZ and pZ→SOZ connections in baseline epileptogenic networks are relevant and, in a two-node model, can be identified with exhaustive BMR rather than specifying models for BMR.

Both of the groupwise covariates evaluated were revealing. Engel class Ia versus Ib outcomes showed greater inhibitory signal from pZ→SOZ. Assuming that inhibitory pZ→SOZ signal is suppressive, then this effect may indicate that not only is SOZ ablation needed for good outcome, but also a stronger pre-operative pZ inhibition encourages possible residual SOZ activity to remain below the seizure threshold. Notably, the potential of measuring a precision network-dependent factor of pZ inhibition to predict seizure outcome from rs-fMRI alone is novel. Further work is indicated, as despite very strong posterior evidence, LOO cross-validation results were not significant. Likely, broader Engel class I-IV outcomes would clarify if the pZ inhibition strength is associated with surgical outcomes.

Absolute signal strength emitted from SOZ and pZ was smaller as age increased. However, these associations did not survive LOO cross-validation, limiting current predicative power of directional signal strength. Decreasing excitation from SOZ→pZ with age is novel, though counterintuitive, given that strength of aberrant signal is expected to increase with disease length. However, recent DCM of epileptiform propagation pathways of Lennox-Gastaut Syndrome also showed reduction with age (Warren et al., 2019). This being noted, the strength difference did not have an effect on SOZ detection yield in this group. Further, the strength difference over age is relatively small, thus not expected to have an effect on SOZ detection in other DRE populations.

Engel outcomes and pre-operative seizure frequency and ROI size were unrelated. Further, age, ROI size, and pre-operative seizure frequency were unrelated to each other. The lack of relationship between surgical outcome and pre-operative seizure frequency and ROI, as well as between variables, make it unlikely that these factors biased Directionality. RSEC strength correlated with Engel outcome, but not seizure frequency nor reduction. Thus, severity of epilepsy also does not appear to influence the applicability of Directionality.

5. Conclusions

This study demonstrates the high accuracy (93%) and sensitivity (96%) of Directionality to identify the origin of excitatory and inhibitory signal between SOZ and the hypothesized pZ in children with HH-DRE. Directionality can be performed from 20 minutes of rs-fMRI data and does not require the additional cumbersome or expensive procedures, personnel, or equipment needed for EEG-fMRI. It also benefits from a lack of need to capture epileptiform activity and ultimately leads to faster results due to greater simplicity of solution – which may also lead to a higher clinical yield. This method validation study in a homogenous population with known SOZ location and surgical outcomes may be helpful in narrowing the SOZ in regions of suspicion for epileptogenicity and may ultimately decrease dependency on iEEG in other DRE populations and other neurological disorders.

Data Availability Statement

The data made available under institutional IRB approval that supports the findings of this study are available in the supplementary material of this article.

CRediT Author Statement

Varina Boerwinkle: Conceptualization, Methodology, Validation, Formal Analysis, Investigation, Resources, Writing – Original Draft, Writing – Review & Editing, Visualization, Supervision, Project Administration; Bethany Sussman: Conceptualization, Methodology, Software, Validation, Formal Analysis, Investigation, Writing – Original Draft, Writing – Review & Editing, Visualization; Sarah Wyckoff: Conceptualization, Investigation, Formal Analysis, Data Curation, Writing – Review & Editing; Iliana Manjón: Visualization, Writing – Review & Editing; Justin Fine: Software, Writing – Review & Editing; P. David Adelson: Writing – Review & Editing

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Declaration of Competing Interest

None of the authors has any competing interests to declare.