Detrended fluctuation analysis in the presurgical evaluation of parietal lobe epilepsy patients

Objective: To examine the usability of long-range temporal correlations (LRTCs) in non-invasive localization of the epileptogenic zone (EZ) in refractory parietal lobe epilepsy (RPLE) patients. Methods: We analyzed 10 RPLE patients who had presurgical MEG and underwent epilepsy surgery. We quantiﬁed LRTCs with detrended ﬂuctuation analysis (DFA) at four frequency bands for 200 cortical regions estimated using individual source models. We correlated individually the DFA maps to the distance from the resection area and from cortical locations of interictal epileptiform discharges (IEDs). Additionally, three clinical experts inspected the DFA maps to visually assess the most likely EZ locations. Results: The DFA maps correlated with the distance to resection area in patients with type II focal cortical dysplasia (FCD) ( p < 0 : 05), but not in other etiologies. Similarly, the DFA maps correlated with the IED locations only in the FCD II patients. Visual analysis of the DFA maps showed high interobserver agreement and accuracy in FCD patients in assigning the affected hemisphere and lobe. Conclusions: Aberrant LRTCs correlate with the resection areas and IED locations. Signiﬁcance: This methodological pilot study demonstrates the feasibility of approximating cortical LRTCs from MEG that may aid in the EZ localization and provide new non-invasive insight into the presurgical evaluation of epilepsy. (cid:1) 2021 International Federation of Clinical Neurophysiology. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
More than 30% of all patients with epilepsy continue to have seizures despite medication (Del Felice et al., 2010).Epilepsy surgery offers an important treatment option for these patients (de Tisi et al., 2011).However, the success of surgery is heavily dependent on the accuracy of presurgical localization of the epileptogenic zone (EZ).The cornerstones of presurgical evaluation are seizure semiology, magnetic resonance imaging (MRI), positron emission tomography (PET), and scalp electroencephalography (EEG) that may be complemented with magnetoencephalography (MEG).The final delineation of the EZ is in many cases performed with intracranial EEG, such as stereotactic EEG (SEEG) (Vakharia et al., 2018), which requires a reasonably accurate hypothesis of the EZ to guide electrode implantation (Kovac et al., 2017).Non-invasive localization of the EZ prior to invasive studies is a genuine challenge for the conventional visual EEG reading.It is especially challenging in extratemporal epilepsies, such as parietal lobe epilepsy, because these cases often exhibit barely recognizable EEG seizure activity or a very fast spread of ictal activity to other cortical areas (Taylor, 2003;Beleza and Pinho, 2011;Ristic et al., 2012;Salanova, 2012;Francione et al., 2015).
In order to support the identification of EZ with interictal SEEG and EEG/MEG data, novel computational tools have been used to quantify various features in the spontaneous EEG activity, such as high frequency oscillations (Jacobs et al., 2018), functional connectivity (Staljanssens et al., 2017), or metrics of scale-free dynamics (Ramon et al., 2008).The scale-free brain dynamics are characterized by power-law long-range temporal correlations (LRTCs), and they are ubiquitous in neuronal activity (Linkenkaer-Hansen et al., 2001).They are often estimated by detrended fluctuation analysis (DFA) (Peng et al., 1994).The scaling exponents of the LRTCs are frequency-band specific, individually characteristic measures that are modulated by stimuli and tasks (Linkenkaer-Hansen et al., 2004;Nikulin and Brismar, 2005;Palva et al., 2013).The LRTCs are an emergent feature of a system that operates near a critical state (Linkenkaer-Hansen et al., 2001;Hardstone et al., 2012;Chialvo, 2010).Indeed, the presence of LRTCs in brain dynamics (Linkenkaer-Hansen et al., 2001), as well as the presence of neuronal avalanches (Beggs and Plenz, 2003) and the correlations among avalanches and LRTCs (Palva et al., 2013;Zhigalov et al., 2015) suggest that the brain operates near a critical state (Chialvo, 2010;Beggs and Plenz, 2003;Plenz and Thiagarajan, 2007;Werner, 2010).The dynamic state is controlled at least by a balanced neuronal excitation and inhibition (E/I) (Beggs and Plenz, 2003;Poil et al., 2012;Meisel et al., 2015;Meisel et al., 2016).Aberrant LRTCs have been reported in many brain disorders, including depression (Lee et al., 2007), schizophrenia (Nikulin et al., 2012), Alzheimer's disease (Stam et al., 2005;Montez et al., 2009), Parkinson's disease (Hohlefeld et al., 2012), and epilepsy (Monto et al., 2007;Ramon and Holmes, 2013), which among other lines of evidence, has lead to the hypothesis that the imbalance of E/I plays a role in the pathophysiological mechanisms of a variety of brain disorders.A shift in E/I balance towards excitation can be seen as an increase in the LRTCs (Poil et al., 2012) until the brain surpasses from the subcritical state (Priesemann et al., 2014) to the supercritical state after which the LRTCs begin to decrease.Excessive excitability is also seen in epilepsy, and indeed, previous interictal EEG and MEG studies have shown an increase in the LRTCs in the vicinity of epileptic brain areas (Parish et al., 2004;Monto et al., 2007;Witton et al., 2019).Parish et al. (2004) showed with SEEG that in unilateral mesial temporal lobe epilepsy the affected hippocampus exhibits higher LRTCs than the contralateral hippocampus.Similarly, an intracranial grid EEG study by Monto et al. (2007) found higher LRTCs in the vicinity of the neocortical EZ compared to the rest of the measured cortex.More recently, a case study by Witton et al. (2019) showed an increase in the LRTCs in epileptogenic brain areas with MEG.The prospect of non-invasive application of the LRTCs to localize the EZ is an appealing novel method to supplement conventional methods in the presurgical evaluation of epilepsy.However, non-invasive studies are still scarce and it is unclear how the LRTCs correlate with more conventional noninvasive methods for EZ localization.
Here, we aimed to determine whether LRTCs estimated by DFA from presurgical interictal MEG recordings correlate with conventionally defined EZ in patients with refractory parietal lobe epilepsy.We also assessed whether DFA scaling exponents correlate with the conventionally identified locations of interictal epileptiform discharges (IED).Additionally, we assessed whether the visual assessment of the spatial maps of the DFA scaling exponents could be used to identify the EZ in a manner that is analogous to conventional radiological assessment.We aimed to emulate a situation where the DFA scaling exponent maps were assessed as any other radiological imaging modality images as per the clinical pipeline.

Patients
Helsinki University Hospital epilepsy surgery register was searched for patients that had undergone resective parietal lobe epilepsy surgery during a 25-year time period between 1991-2016.The inclusion criteria for this study were that the patient had: (1) undergone resective parietal lobe epilepsy surgery; (2) a minimum of two years of post-operative follow-up; (3) preoperative MEG recordings with a minimum of 10 min of continuous data without ictal events; (4) no previous cortical resection or large cortical malformations; (5) availability of pre-and postoperative T 1 -weighted MR images.
In total 10 patients fulfilled these criteria.25 patients were identified of whom 18 had preoperative MEG recordings.The MEG data of one of these patients was irretrievable.Out of the 17 patients left, 14 had at least 10 min of continuous MEG data.Four patients were additionally excluded due to a previous cortical resection, a large cortical malformation that extended over most of the affected parietal lobe complicating segmentation, absence of post-operative MRIs, or excessive edema in the MRI that precluded accurate delineation of resection area.
All patients had multiple MEG measurements with different conditions performed as per clinical protocol.In these cases, all measurements that complied with the inclusion criteria were included.In total 16 MEG measurements of 10 patients were included (Table 1).
We grouped patients according to etiology into three categories, namely "FCD I", "FCD II", and "Other".These groups contained 3, 4, and 3 patients respectively (Table 1).The Other group included patients with perinatal infarction, intraventricular hemorrhage, and one patient with simultaneous perinatal infarction and FCD type III.
The study was approved by the HUS Medical Imaging Center, University of Helsinki and Helsinki University Hospital.All data were collected retrospectively from the patient records and from the studies which had been carried out as part of the routine preoperative evaluation.

Measurements
Electrophysiological recordings were performed with 306channel Elekta Vectorview TM (MEGIN (Elekta Oy), Helsinki, Finland) consisting of 204 planar gradiometers and 102 magnetometers at 600 Hz sampling rate.An online digital band-pass filter of 0:1 to 200 Hz was applied in every recording.Head position indicator (HPI) coils were used to determine the head position within the MEG gantry.The HPI coils were referenced to three anatomical fiducial points (left and right preauricular points, and nasion) before the MEG measurements.Seven patients had continuous head position tracking during the measurements.

Data preprocessing
The raw MEG data were processed with tSSS using MaxFilter software (MEGIN) (Taulu et al., 2005;Taulu and Simola, 2006).We inspected the recordings visually to exclude bad channels and ictal events.Ocular and heartbeat artifacts were detected with electrooculogram and electrocardiogram respectively when available.We isolated these components with signal-space projection as implemented in Brainstorm, and they were subsequently removed.
We band-pass filtered the data with finite impulse response (FIR) filter to four bandwidths of 1:00; 6:25 ½ Hz, 4:00; 16:00 ½ Hz, 11:00; 44:00 ½ Hz, and 3:00; 48:00 ½ Hz, which are here referred to as delta, alpha, beta, and broadband (BB) bands respectively.The FIR filter was based on a Kaiser window design as implemented in Brainstorm.The ripple in pass and attenuation in stopband were set to 10 À3 and 60 dB respectively.We excluded the band-pass transients at the beginning and the end of the data from further analysis.

Source projection and surface parcellation
Pre-operative anatomical MR images were processed with the FreeSurfer software suite (https://surfer.nmr.mgh.harvard.edu/).For individual MR protocols see Supplemental Information.The processing pipeline included volumetric segmentation, surface reconstruction, and cortical parcellation and labeling with Schaefer 2018 atlas (Schaefer et al., 2018) with the resolution of 200 parcels.
We performed the MEG-MRI coregistration and all subsequent analysis with Brainstorm (Tadel et al., 2011) MATLAB toolbox (MathWorks, Natick, MA, USA).The three fiducial points were automatically placed in Brainstorm on the individual MRIs and MEG was coregistered by corresponding points.Further refining was performed using digitized head surface points.One patient (patient 02) did not have a digitized head surface, and in their case, only fiducial points were used.
Brainstorm was used to compute overlapping spheres forward models and cortically constrained source models.We performed cortically constrained source estimation with the sLORETA method (Pascual-Marqui, 2002) to estimate time series for 15 000 source vertices with an average vertex resolution of 4:4 mm.Due to the lack of empty room measurements, identity matrices were used as noise covariance matrices.The source vertices were then collapsed into parcel times series by calculating the average of the sign-flipped time series of the vertices within the parcel.The sign flipping depended on the source orientation so that the sign of the time series of a source was flipped if the orientation of the source was against the dominant orientation.We defined the dominant orientation of a parcel by singular value decomposition of all source dipoles within the parcel.

Detrended fluctuation analysis
Detrended fluctuation analysis was used to estimate the LRTC scaling exponents (Peng et al., 1994;Linkenkaer-Hansen et al., 2001).DFA implementation was adapted from the implementation published earlier by Hardstone et al. (2012).First, we calculated the amplitude envelope A t ð Þ via the Hilbert transform.Then, we b extent of the resection limited by the proximity of the primary sensorimotor cortex.c based on the effect on cognition; no significant reduction in seizures.
calculated the cumulative sum of the amplitude envelope to get the signal profile: We divided the signal profile (Y t ð Þ) into W windows of size s with 50 % overlap.A linear trend within each window was calculated using a least-squares fit and subsequently removed from the signal profile: We then computed the fluctuation F s ð Þ by taking the standard deviation of the detrended signal profile Y dtr t; s ð Þ for each window: where N s ð Þ is the number of samples within the time window of size s.The median of the fluctuation F s ð Þ h iwas calculated and this process was repeated for 20 logarithmically spaced time windows between 10 and 60 s.
The rationale for using median instead of mean as was used in Hardstone et al. (2012) was to exclude the possible effects of largeamplitude artifacts that may bias the mean.
We then plotted the median fluctuation values against the time windows on logarithmic axes (Fig. 1B).These quantities were linearly correlated and the slope was defined as the DFA exponent a: We use a as an estimate for the LRTC scaling exponent.
We calculated DFA exponents for each parcel time series of each frequency band of each measurement.With 16 measurements all band-pass filtered into 4 frequency bands, this resulted in 64 DFA exponent distributions.The DFA exponents were then mapped onto the pial surface of each patient to form DFA maps for visualization.

Interictal epileptiform discharge analysis
IED detection and localization were performed separately in an independent study by Wilenius et al. (2020).In brief, thirty minutes of interictal MEG data were selected from the measurements of each patient.The MEG epochs used for the IED detection were generally longer and only partly overlapping the epochs used for the DFA.IEDs were identified visually from the MEG data and classified to spike types by an experienced clinical neurophysiologist (JW).Further IED identification based on the spike types were performed with a pattern search function in BESA Research 6.1 (BESA GmbH, Gräfelfing, Germany).The IEDs of each patient were averaged based on the spike type.Equivalent current dipole with a single-shell spherical conductor model was performed using the FieldTrip MATLAB toolbox (Oostenveld et al., 2011).The time interval for dipole fitting was selected to start at the onset of the discharge and to end at the earliest peak of the spike.

Visual assessment of DFA maps by clinical experts
Three experienced clinical neurophysiologists (LL, MP, and SV) reviewed the DFA exponent distributions projected on individual 3D cortical images.The DFA maps were shown to them in the same session without any other clinical information.Everyone gave their assessment independently and blinded from each other's assessment.With the assumption that the DFA exponent is increased near the EZ, they were tasked to visually assess the most likely hemisphere and lobe for the EZ, based solely on the spatial distribution of the DFA exponents.Then, the experts were tasked to assess the most likely hemisphere and sub-lobe (superior or inferior) given that the lobe is parietal in these patients.In ambiguous cases, the experts had an option to answer "I do not know", which was treated as missing data.

Data analysis and statistics
We assessed the relationship between the DFA exponents and distance to the resection area and the IED locations by linear regression and correlation coefficients.Since we sought to examine the intrapatient correlations between the DFA exponent distribution and the resection area while maintaining ease of comparison between patients, we rescaled the DFA exponents linearly between 0 and 1 for each patient and frequency band.We defined the normalized DFA exponent a 0 as: where A is the set of DFA exponents of one measurement at one frequency band, and a is a member of A.
We defined distance as the shortest distance between the centroids of adjacent parcels (Dijkstra, 1959).The results were grouped by the three etiologies and by frequency bands.We calculated the coefficient of linear regression and the Pearson correlation coefficient for the first 10 cm from the most central parcel of the resection area and the IED location to the centroid of every other parcel.Ideally, the change in the DFA exponents would be very focal in the vicinity of the resection area.This implies that including the whole hemisphere in larger brains would bias the linear regression towards the baseline levels.To avoid this bias, linear regression was estimated only up to a standard distance from the parcel of interest; it was aimed to cover the largest distance from the resection center to the furthest adjacent non-resected parcel in all patients (5:8 cm), while also being less than the smallest distance from the resection center to the furthest parcel over all patients (17:5 cm).A compromise between these needs was chosen at 10 cm.
The results of the visual assessment were likewise grouped by etiology.We calculated accuracy relative to the resection area and respective interobserver agreement (IOA) for each group and each frequency band separately.The IOA was quantified by using the Krippendorff's alpha (K a ) (Krippendorff, 1970;Krippendorff et al., 2004).K a measures observed disagreement relative to the disagreement expected by chance.It was used instead of Cohen's kappa due to K a allowing more than two observers and being able to deal with missing data.
In order to assess the collinearity of the fluctuation values against time windows and the deviation from the power law, we calculated Pearson correlation coefficient and a parabolic index, b, for each a (Monto et al., 2007).We defined the parabolic index as: where E 1 and E 2 are the mean-squared errors (MSE) of the 1st-and 2nd-order least-squares fits of log F s ð Þ h i against log s ð Þ respectively (Fig. 2A).The values of b range from 0 to 1.If the data follows a linear trend, i.e. the fluctuations obey the power law, the 1st-order MSE is low and b tends to 0. High b-values could indicate that the data does not conform to power law.
In order to inherently and rigorously account for multiple comparisons, we generated surrogate data by randomly permuting the locations of the resection centers and the IED location of a randomly chosen patient, and performing the aforementioned analyses on the surrogate data.For each set of surrogate data, we took the greatest magnitude statistic between the four bandwidths as the maximal statistic.We repeated this until all possible permutations or a sufficient sample of different permutations were included to generate maximal statistic distribution against which the experimental statistics could be compared to (Nichols and Holmes, 2002).For further details on the maximal statistic approach, see Supplemental Information.The 95%-confidence

Table 2
Results of the visual inspection by three experts.The experts were tasked to independently mark the hemisphere and the lobe that visually appeared to have the highest detrended fluctuation analysis (DFA) exponents.Then, given the correct lobe but not hemisphere, they were to mark the hemisphere and sub-lobe (superior or inferior part of the lobe) based on the same criterion.The uncertainty of accuracy is the 95%-confidence interval.The interobserver agreement is quantified by Krippendorff's alpha.In bold are those accuracy-values that statistically significant (family-wise error rate -adjusted two-tailed p < 0:05).The asterisk signifies those values where the significance level is at the minimum.Abbreviations: FCD, focal cortical dysplasia.intervals for accuracy of visual inspection analysis were calculated by bootstrapping (Efron, 1979).The usage of surrogate data enables estimating the chance level without assumptions on the distribution of data, and using the maximal statistic approach controls the family-wise error rate (FWER) (Nichols and Holmes, 2002).

Visual assessment of DFA maps
The visual assessment of the DFA maps showed typically one or a few clusters of higher DFA values.Hence, akin to the clinical visual assessment of any processed neuroimaging data, we sought out to first determine whether a visual assessment of the DFA maps could provide a reliable prediction of the EZ.To this end, three experts assessed the DFA maps from all frequencies, while being blinded to each other as well as to the patient information.We estimated the accuracies, related p-values, and interobserver agreements for all frequency bands (Table 2).

Hemisphere and lobe
The correct hemisphere was identified with an accuracy of 0:69 AE 0:18 (p < 0:05) and an interobserver agreement (IOA) of 0:82 when averaging over all frequency bands and etiologies.Etiology wise, significantly above chance level lateralization accuracy and high interobserver agreement was achieved only in the FCD groups, whereas in the other etiologies the accuracy was at chance Fig. 4. Relation between detrended fluctuation analysis (DFA) exponent and distance from interictal epileptiform discharge (IED) location.The linear regression was performed with normalized DFA exponents and non-normalized distance.Patients with in "FCD I" and "Other" categories had multiple different IED locations.In these cases the IEDs with overall best agreement with DFA exponent is shown.The leftmost y-axis represents the mean maximum, mean, and minimum non-normalized DFA exponent values for each bandwidth.Binning was performed after the regression analysis.Each bin is 1 cm wide.The whiskers correspond to 1:5 times the interquartile range.The significance is tested against the resection center being in any other parcel.The significance level is represented by family-wise error rate -adjusted two-tailed p-value.Abbreviations: FCD, focal cortical dysplasia.level (average accuracy: 0:51 AE 0:27).The two FCD groups were not significantly different in lateralization accuracy (p > 0:05).
The correct hemisphere and the correct lobe were identified with an accuracy of 0:29 AE 0:19 (p < 0:05) and IOA of 0:59 over all patients.The change-level accuracy was 0:125.The accuracy was better in the FCD groups than in the other etiologies (0:41 AE 0:23 and 0:11 AE 0:16 respectively) and the two FCD groups were again not significantly different (p > 0:05).

Hemisphere and sub-lobe
When lateralizing the known parietal lobe, the lateralization accuracy was slightly poorer than when the lobe was not known (0:61 AE 0:22; p > 0:05).However, this difference was not significant (p > 0:05).High lateralization accuracy was achieved in the combined FCD group but not in the other etiologies (0:79 AE 0:14 and 0:39 AE 0:32 respectively).The two FCD groups were not significantly different in accuracy (p > 0:05) when compared to each other.
Visually identifying both the correct hemisphere and the correct sub-lobe (inferior or superior) of the known lobe was feasible only in the delta band (Table 2) when all etiologies were taken into account.Out of the three etiology groups, the FCD II group was the only group where both the correct hemisphere and the correct sublobar region of the known lobe ware consistently identified (0:73 AE 0:24; p < 0:05; IOA ¼ 0:80, chance-level accuracy: 0:25).

Distribution of DFA exponents against resection area and IED locations
The DFA exponents had a significant (p < 0:05) negative linear correlation with the distance to the resection area in all frequency bands in the FCD II patient group (Fig. 3).A similar negative correlation was observed in the FCD I etiology group, although the effect was less pronounced and statistically significant only at the beta band.The group with the other etiologies showed no correlation between the DFA exponents and the distance to the resection area.The magnitude of slope and correlation coefficients were the highest in delta and broadband bands over all patients.
All patients with FCD II had one IED location per patient and their IED location correlated strongly with the DFA exponent distribution.Concurrently, the DFA exponents also correlated strongly with their IED locations in all but beta frequency band (Figs. 4  and 5).
Among the 3 FCD I patients, 7 distinct IED locations were observed of which 5 localized on the operated hemisphere.We observed no correlation between the DFA exponents and the IED locations in FCD I patients (Figs. 4 and 5).
Similarly, among the 3 patients with etiologies other than FCD, we observed 8 distinct IED locations.3 of these localized on the operated hemisphere.We observed no correlation between the DFA exponents and the IED locations in this etiology group.

Discussion
Our work shows that the LRTCs are increased near the epileptogenic zone defined by the surgical resection area.This relationship was observed both with correlation analyses and in a clinical visual inspection of the cortical DFA exponent maps.The findings are in line with previously reported results (Parish et al., 2004;Monto et al., 2007;Witton et al., 2019).However, here we extend the prior results by presenting the analysis using non-invasive recordings in a clinically difficult-to-assess group of extratemporal epilepsies.Our data also suggest that relationship between the cortical LRTCs and the EZ may be dependent on etiology, which may offer new insights into the understanding of the neural network mechanisms associated with different etiological subgroups.
We found the highest and the most focally increased DFA exponents in patients who had type II FCD, most of whom also had Engel I outcome in the two-year post-surgery follow-up.Additionally, three out of four of the FCD II patients were MRI-negative (Table 1).The DFA exponent was also observed to be increased near the resection area in patients with type I FCD, although the effect was not as pronounced as with FCD II patients.All FCD I patients also had worse outcomes (III and IV) than the FCD II patients.In patients with perinatal infarctions or double etiologies, no correlation between the DFA exponents and the resection area was observed.These findings were corroborated by visual inspection by expert clinical neurophysiologists.The hemisphere was accurately and with high interobserver agreement recognized from the DFA maps only in the FCD patients (both type I and type II).The lobe or sub-lobe were reliably identified only in FCD type II patients.With respect to the frequency band, the alpha and the delta bands showed the highest congruence to the resection area, whereas the beta band the worst.In previously reported intracranial EEG results (Monto et al., 2007) the beta band was observed to be the most sensitive to epileptic focus.However, detecting high-frequency activity is more difficult with MEG than with intra-cranial methods due to poor signal-to-noise ratio (Dalal et al., 2009;Jerbi et al., 2009).Yet, we did observe statistically significant correlation between the beta band DFA exponents and the resection area in FCD patients, however, this did not materialize in the visual assessment.
There are some limitations that affect the generalizability of the present results.The cohort was relatively small and retrospective, and we only included patients with parietal lobe epilepsy that had undergone epilepsy surgery.A technical limitation in our data set is the heterogeneity.The cohort included different etiologies with a limited number of subjects in each category.In addition, patients' state could not be fully controlled in this retrospectively collected cohort.Data segments included tactile or electric somatosensory stimulations, hyperventilation, or sleep and awake resting-state measurements.Linkenkaer-Hansen et al. (2004) showed that external somatosensory stimulation decreases the magnitude of the DFA exponent, particularly in the somatosensory cortex.On the other hand behavioral state (sleep vs. awake) has been shown to globally change the magnitude of the DFA exponent (Kim et al., 2009;Lee et al., 2002).Anti-epileptic drug (AED) use was also not controlled for (Meisel et al., 2015;Meisel et al., 2016).AEDs are known to alter cortical dynamics and thus affect the DFA exponent (Monto et al., 2007;Poil et al., 2011).Notably, these prior findings report global effects on DFA by situational factors, hence they do not directly challenge our present work showing spatially selective DFA distributions.While these data provide proof of concept for the potential of DFA to localize epileptogenic zone, these data do not enable the estimation of the actual added value to the clinical planning process of epilepsy surgery.Further studies are needed with unselected and more variable patient population, including such patients that were not eventually selected for surgery.
The apparent correlation with etiology is an interesting observation and solicits a closer investigation.To our knowledge, there have not been previous reports on the etiological dependence of the DFA exponents related to epilepsy.FCD II is known to exhibit distinct interictal activity compared to other etiologies within the cortical lesions as measured with SEEG (Tassi et al., 2012;Di Giacomo et al., 2019).The EZ in the FCD II also tends to be more focal whereas in the FCD I and in other etiologies it often involve multiple lobes (Xue et al., 2016;Isler et al., 2017).This corroborates our observations that the DFA exponent distributions were more focal in the FCD II patients than in the other patients.
While promising, this is primarily a methodological pilot study and there are several checkpoints that DFA-based assessment has to pass before wider clinical adoption.At first, the specific use case for DFA among the existing methodological repertoire in presurgical planning needs to be defined (Mouthaan et al., 2016).An obvious clinical need for such a method is in planning SEEG implantation or subdural grid electrode placement.This requires sufficient, but not necessarily perfect, accuracy in non-invasive seizure onset zone localization.Our DFA pipeline is transparent and automatizable.It may also be applied to high-density EEG recordings, preferably using source-level signals (Lai et al., 2018;Tokariev et al., 2016).Our work showed quantifiable distribution patterns speaking for an automated localization.Additionally, DFA may be used in a more classical manner, akin to any common radiological assessment: the DFA exponent maps can be assessed visually, as was done in our present work.This could lift some of the concerns that clinicians might have with the use of fully automated analysis methods.Finally, the critical question before wider clinical adoption is whether clinical experts will see a perceived added value from DFA over other, already existing non-invasive methods.

Conclusions
DFA exponent as a measure of LRTC is a promising tool for presurgical evaluation of epilepsy.It may aid in localizing the epileptogenic zone from spontaneous EEG, thus providing a complementary tool to IED localization in interictal recordings.We found the LRTC in the epileptogenic zone to correlate with the epilepsy etiology.As this was primarily a methodological pilot study, further exploration of this observation is warranted and might offer additional insight into the pathological dynamics in epilepsy and thus provide new information of epilepsy, in addition to potentially improving patient care.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Overview of data acquisition and analysis.(A) MEG measurements and MRIs, and the relevant clinical information of 10 patients with drug-resistant parietal lobe epilepsy were acquired as part of clinical protocol.(B) The MRIs were used to create individual head models, which were used in the calculate of cortical source signals from bandpass filtered MEG.The cortical surface was parcellated into 200 parcels and the source signals were collapsed into these parcels.detrended fluctuation analysis (DFA) was performed for each parcel signal resulting in cortical DFA exponent (a) distribution.(C) These DFA exponent distributions were visualized on the cortical models of the patients for quantitative and visual analyses.In the quantitative analysis the distance from the resection area, the distance from interictal epileptiform discharge (IED) locations, and the magnitude of the DFA exponents were correlated and cross-correlated.In the visual analysis three experienced clinical neurophysiologists assessed the 3D cortical DFA distribution visualizations and determined the most likely location of seizure onset zone based on the DFA exponent distributions.We then grouped the results of these analyses by the etiology of the patients.

Fig. 2 .
Fig. 2. (A) Example figures demonstrating the collinearity of median fluctuation (F) and time window (s) at different parabolic indices (b).b-value was used to assess the deviation from power law behavior (see Methods).Each example represents data from one parcel at one frequency band of one measurement.The three b-values correspond to 50%; 75%, and 95% cut-points of all b-values.Low b-value indicates power law behavior.(B) The distribution of b-values at different frequency bands.(C) The relationship between Pearson correlation coefficient (r) and b-value.(D) The distribution of the detrended fluctuation analysis (DFA) exponents (a) in different MEG measurements.The gray zone indicates the range of DFA a values that indicate long-range temporal correlations.The horizontal axis shows the measurement number and the corresponding patient number.Some patients had multiple MEG measurements performed as part of the clinical routine.The Pearson correlation coefficients (r) above the repeated measurements indicate the repeatability of each measurement pair.

Fig. 3 .
Fig. 3. Relation between detrended fluctuation analysis (DFA) exponent and distance from resection center at different bandwidths and different etiology groups.The linear regression was performed with normalized DFA exponents and non-normalized distance.The leftmost y-axis represents the mean maximum, mean, and minimum nonnormalized DFA exponent values for each bandwidth.Binning was performed after the regression analysis.Each bin is 1 cm wide.The whiskers correspond to 1:5 times the interquartile range.The significance is tested against the resection center being in any other parcel.The significance level is represented by family-wise error rate -adjusted two-tailed p-value.Abbreviations: FCD, focal cortical dysplasia.

Fig. 5 .
Fig. 5. (A) Pearson correlation coefficients (r) and partial correlation coefficients between three variables, namely detrended fluctuation analysis (DFA) exponent (a), distance to the resection center (Ddr), and distance to interictal epileptiform discharge (IED) location (Dds).In each patient, the IED location closest to the resection area center was considered.The correlations are mean correlations over all frequency bandwidths and the uncertainty is the corresponding standard deviation.Values in bold red are statistically significant (two-tailed p-value < 0:05).(B) Two examples (Patient 04 Spontaneous and patient 05 MN SEF) of DFA exponent distributions at broadband bandwidth.The green outlines show the extent of the surgical resection.Abbreviations: FCD, focal cortical dysplasia.
a MEG study during focal status epilepticus; FPHT and VPA given as loading doses.