## Abstract

Epilepsy is a neurological disorder characterized by the sudden occurrence of unprovoked seizures. There is extensive evidence of significantly altered brain connectivity during seizure periods in the human brain. Research on analyzing human brain functional connectivity during epileptic seizures has been limited predominantly to the use of the correlation method. However, spurious connectivity can be measured between two brain regions without having direct connection or interaction between them. Correlations can be due to the apparent interactions of the two brain regions resulting from common input from a third region, which may or may not be observed. Hence, researchers have recently proposed a sparse-plus-latent-regularized precision matrix (SLRPM) when there are unobserved or latent regions interacting with the observed regions. The SLRPM method yields partial correlations of the conditional statistics of the observed regions given the latent regions, thus identifying observed regions that are conditionally independent of both the observed and latent regions. We evaluate the performance of the methods using a spring-mass artificial network and assuming that some nodes cannot be observed, thus constituting the latent variables in the example. Several cases have been considered, including both sparse and dense connections, short-range and long-range connections, and a varying number of latent variables. The SLRPM method is then applied to estimate brain connectivity during epileptic seizures from human ECoG recordings. Seventy-four clinical seizures from five patients, all having complex partial epilepsy, were analyzed using SLRPM, and brain connectivity was quantified using modularity index, clustering coefficient, and eigenvector centrality. Furthermore, using a measure of latent inputs estimated by the SLRPM method, it was possible to automatically detect 72 of the 74 seizures with four false positives and find six seizures that were not marked manually.

## 1 Introduction

### 1.1 The Caveat in the Correlation, Precision Matrix, and Sparse-Regularized Precision Matrix Methods

The correlation method is widely used for estimating brain functional connectivity (Anand et al., 2005; Biswal, Yetkin, Haughton, & Hyde, 1995; Bruno et al., 2017; Chu et al., 2015; Lynall et al., 2010; Rubinov & Sporns, 2010; Saggar et al., 2015; Shafi, Westover, Oberman, Cash, & Pascual-Leone, 2014; Siegle, Thompson, Carter, Steinhauer, & Thase, 2007; Smith et al., 2011; Uddin et al., 2009; Vertes et al., 2012; Zhou, Thompson, & Siegle, 2009). This method has been used to analyze large-scale resting-state brain networks (Biswal et al., 1995; Cordes et al., 2000; Fox, Corbetta, Snyder, Vincent, & Raichle, 2006; Fox et al., 2005; Greicius, Krasnow, Reiss, & Menon, 2003; Hampson, Peterson, Skudlarski, Gatenby, & Gore, 2002; Uddin et al., 2009; Xiong, Parsons, Gao, & Fox, 1999) in healthy subjects, as well as altered brain functional connectivity in pathology (Anand et al., 2005; Supekar et al., 2008) in the human brain. However, connectivity estimation using the correlation method can be misleading or incorrect since brain regions might show high correlation due to a common input, which may or may not be measured or observed, and not due to strong physical connections between them (Wang, Kang, Kemmer, & Guo, 2016). In a recent functional magnetic resonance imaging (fMRI) study from 210 healthy young adults in the Human Connectome Project (HCP), researchers have been able to identify many new areas in the human cortex using a machine learning classifier (Glasser et al., 2016). Pairs of these areas with a high degree of functional connectivity also received common input from other areas of the human brain. Honey et al. (2009) measured the functional connectivity from fMRI in five individuals using the correlation method, and several brain regions showed a high degree of functional connectivity without having direct connections. This provides further evidence that brain regions can have high correlations due to common inputs rather than direct interactions.

Due to this major drawback of the correlation method, researchers have suggested using partial correlations to measure strong, direct interactions between pairs of brain regions while simultaneously removing the influence of the rest of the brain regions (Das et al., 2017; Dempster, 1972; Lauritzen, 1996; Whittaker, 1990), assuming all the regions can be measured. Partial correlations help identify pairwise brain regions that are conditionally independent given all the other brain regions. When the output of the brain regions follows a multivariate gaussian distribution, the inverse covariance matrix (ICM) or precision matrix (PM) can be used to calculate these partial correlations (Dempster, 1972). A value of zero or very close to zero in the PM will indicate that the two brain regions are conditionally independent given the rest of the brain regions. Researchers (Marrelec et al., 2007; Salvador et al., 2005) have previously used the PM to estimate the human brain functional connectivity from fMRI.

Hsieh et al. (2011, 2013) have applied the SRPM method to an fMRI data set collected from the human brain. After applying modularity-based clustering (Blondel, Guillaume, Lambiotte, & Lefebvre, 2008; Brandes et al., 2008; Newman & Girvan, 2004; Newman, 2006; Reichardt & Bornholdt, 2006; Ronhovde & Nussinov, 2009; Sporns, 2010; Sun, Danila, Josić, & Bassler, 2009) to the regularized precision matrix, they found strong, functionally connected regions in gray matter regions in the human brain, along with identification of resting-state networks such as default mode and sensorimotor networks. The modules detected by the SRPM method were similar to those identified using independent components analysis (ICA) on the same data set. Ryali, Chen, Supekar, and Menon (2012) applied the SRPM method to resting-state fMRI data and found a modular architecture characterized by strong interhemispheric links, distinct ventral and dorsal stream pathways, and a major hub in the posterior medial cortex. Varoquaux, Gramfort, Jean-Baptiste, and Thirion (2010) analyzed human brain functional connectivity using the SRPM method; after clustering, they found regions corresponding to important brain areas such as the primary visual system (medial visual areas), the dorsal visual pathway, the occipital pole, and the intraparietal areas comprising the default mode network, the fronto-parietal networks, the ventral visual pathway, the lateral visual areas, and the inferior temporal lobe. Monti et al. (2014) analyzed human brain functional connectivity using the SRPM method in healthy patients and found that the right inferior frontal gyrus and the right inferior parietal lobe play a key role in attention and executive function during cognitively demanding tasks and may be fundamental in regulating the balance between other brain regions.

Rosa et al. (2015) also have applied the SRPM method to distinguish patients with major depressive disorder from healthy control subjects by analyzing brain functional connectivity from fMRI data. Allen et al. (2012) analyzed brain functional connectivity using SRPM and clustering algorithms from resting-state data in young adults and found connections between regions in the lateral parietal and cingulate cortex. They also found that the dynamic functional connectivity of the human brain was markedly different from the stationary brain connectivity. The SRPM method has also been used to estimate brain functional connectivity at a neuronal level in the mouse visual cortex in Yatsenko et al. (2015), where the authors argued that this method found more biologically plausible brain networks than the correlation method did. Wang, Kang, Kemmer, and Guo (2016) applied the SRPM method to estimate human brain functional connectivity from fMRI data and compared it with the connectivity estimated by the correlation method. They found that the SRPM method removed considerable between-module connections that were estimated by the correlation method. In particular, $34%$ of the significant connections found in the correlation method became insignificant after calculating the partial correlations from the SRPM method. Recently, Glasser et al. (2016) used fMRI data of the human cerebral cortex from 210 healthy young adults and found 97 new cortical areas, along with 83 previously existing cortical areas, with the SRPM method. They concluded that their study presents a picture of the structural and functional organization of the human cerebral cortex and its variation across individuals and in development, aging, and disease.

Although the application of the SRPM method to identify functionally connected networks in the human brain is promising, one major shortcoming of this method is that it assumes that all brain regions are measurable or observed. This assumption might lead to incorrect estimation of brain connectivity since most brain regions remain unobserved (also known as latent regions) using current technologies, especially in electrocorticography (ECoG) recordings, which are used in our analysis.

### 1.2 Proposed Method: The Sparse-Plus-Latent-Regularized Precision Matrix Method

### 1.3 Prior Research on Human Brain Connectivity Estimation Using the SLRPM Method

Research on the application of the SLRPM method on brain functional connectivity estimation has been limited, and there is as yet no comprehensive analysis on the application of the SLRPM method for human brain connectivity estimation. Yatsenko et al. (2015) have applied the method to infer connectivity in the mouse visual cortex. Their simulation results showed that the SLRPM method has the potential to outperform both the correlation and the SRPM methods. However, they generated surrogate data from the same optimization problem that they solved to estimate the connectivity matrix for the SRPM and SLRPM methods. A less biased analysis would generate data independent of the methods and then validate the methods by estimating the connectivity matrices from these data. This approach was taken in our simulations for validating the methods. From the experimental analysis using 3D random-access laser scanning microscopy of calcium signals in the mouse visual cortex, the authors found more physiologically interpretable functionally connected brain regions by the SLRPM method than the correlation and SRPM methods. Our effort in this letter on the use of the SLRPM method to characterize brain connectivity in humans from ECoG recordings is the first.

### 1.4 First Contribution: Validating SLRPM Using Artificial Networks

In this letter, we generate artificial networks via simulation and demonstrate the potential of the SLRPM method in recovering network connectivity when only part of the network is observed and the rest of the network is latent or unobserved. For comparison, we consider the correlation method, the precision matrix method, and the SRPM method. This is our first contribution. To evaluate the performance of the SLRPM method, we generated synthetic data from a mechanical model that consists of a cascade connection of a number of springs and masses. The purpose here is to demonstrate that the connectivity pattern of the observed spring-mass network can be obtained from the SLRPM associated with the displacements of the observed masses. We considered several cases: sparse and dense connections, short-range and long-range connections, and varying the number of latent variables. The performance of the SLRPM method was superior to those of the other methods.

### 1.5 Second Contribution: Application of SLRPM for Epileptic Seizure Analysis

Epilepsy is a neurological disorder characterized by the sudden occurrence of unprovoked seizures, affecting more than 50 million people worldwide. There is extensive evidence of significantly altered brain connectivity during epilepsy, especially during seizure periods, in the human brain (Diessen, Diederen, Braun, Jansen, & Stam, 2013; Kramer & Cash, 2012; Stam, 2014). However, research on the characterization of human brain functional connectivity during epileptic seizures has predominantly been limited to the use of the correlation method.

Vega-Zelaya, Pastor, de Sola, and Ortega (2015) estimated brain functional connectivity during partial seizures from foramen ovale electrode (FOE) and electroencephalogram (EEG) recordings in 22 temporal lobe epilepsy (TLE) patients using the correlation method. Brain network connectivity was characterized using the average clustering coefficient (ACC) and the modularity index (MI). The ACC increased and the MI decreased after seizure onset. Kramer et al. (2010) analyzed brain functional connectivity in 48 seizures from ECoG recordings in 11 epilepsy patients using the correlation method. They found that the ACC increased after the seizure onset. A research study by Ponten, Bartolomei, and Stam (2007) characterized, via the correlation method, the functional connectivity from EEG recordings in seven patients suffering from mesial TLE (MTLE) and reported an increase in ACC after seizure onset. In another study, Schindler, Bialonski, Horstmann, Elger, and Lehnertz (2008) carried out an analysis of brain activity from EEG recordings in 100 epileptic seizures from 60 patients using the correlation method. They characterized functional connectivity by ACC and found that it increased after seizure onset. Burns et al. (2014) used another measure to characterize brain functional connectivity in terms of eigenvector centrality (EC), from ECoG recordings in 12 patients before and after seizure onset via the correlation method. They found that brain functional connectivity significantly changed before and after seizure onset, and they successfully localized the seizure onset zones by analyzing this change in connectivity (EC) of the ECoG electrodes.

These results of the human brain functional connectivity during epileptic seizures measured using the correlation method might be flawed since existing brain imaging technologies cannot simultaneously record from the entire human brain, and brain regions that are observed can receive common inputs from latent or unobserved brain regions. Furthermore, there can be dopaminergic, norepinephrine, or cholinergic fiber projections to the cortex (Attwell & Iadecola, 2002; Giorge, Pizzanelli, Biagioni, Murri, & Fornia, 2004; Krimer, Muly, Williams, & Goldman-Rakic, 1998; Lado & Moshe, 2008; Sato & Sato, 1992), and these neuromodulatory projections serve as common latent inputs to the cortex and induce spurious correlations among the EEG or ECoG electrodes. But the SLRPM method calculates the precision matrix of the conditional statistics of the observed brain regions given the latent brain regions, and this helps to remove the influence of common observed and latent inputs on the EEG or ECoG recordings. In this way, we may be able to find more accurate brain connectivity from the recordings, which will help us analyze brain activity during epileptic seizures.

In our second contribution, we demonstrate the application of the SLRPM method to estimate brain connectivity during epileptic seizures from human ECoG recordings. These recordings are analyzed using SLRPM and Louvain method for community detection (LMCD) (Newman, 2006), and brain connectivity has been quantified using modularity index (MI), clustering coefficient (CC), and eigenvector centrality (EC).

### 1.6 Organization of the Letter

In section 2, we test our methods using the spring-mass model and demonstrate the superior ability of the SLRPM method for recovering the connectivity of the spring-mass network when part of the network is unobserved. In section 3, we demonstrate the application of the SLRPM method for estimating brain connectivity during seizures from human ECoG recordings. We draw conclusions in section 4.

## 2 The Spring-Mass Model

However, we assume that we can observe only part of the entire network because some of the displacements are not available for measurement. Hence, the only available data are the sample covariance matrix of the displacements, which are available for measurement. From this sample covariance matrix, we wish to recover the connectivity of the observable masses using the correlation method, the precision matrix (PM) method, the SRPM method, and the SLRPM method. We study the recovery performance of the methods by varying the number of observable variables $p$ (hence, varying the number of latent variables $h$). The values of $p$ are taken to be 60, 100, and 150 for this example. This example is shown in Figure 2. The connections in the ground-truth matrix of the observed network are shown in black dots. We choose the $M$ largest elements (in absolute value) in the estimated connection matrix, excluding the diagonal, from each method, where $M$ is the number of connections in the ground-truth matrix, excluding the diagonal, and check whether the estimated connection corresponds to the connection in the ground-truth matrix. For each method, we show a correctly identified connection in black dots and an incorrectly identified connection in red dots in Figure 2. The error percentage for a particular method is the percentage of incorrectly identified connections expressed as a fraction of the total number of connections in the ground-truth matrix.

Observe that the error percentage that we report here is equal to the number of false positives, which is also equal to the number of false negatives in this case since we only select $M$ largest elements. For example, if we assume that there are five ground-truth connections and the number of hits (true positives) by a method is three, then the number of false negatives is two since the number of misses is two and the number of false positives is also two since those two false connections will appear elsewhere. The value of the regularization parameter $\lambda $ for the SRPM method and the values of the regularization parameters $\alpha $ and $\beta $ for the SLRPM method are also shown in Figure 2. These values were fine-tuned for this example. From this figure, we observe that the PM method, which directly inverts the sample covariance matrix, results in extremely poor performance for the three cases considered. The SLRPM method has comparable or better performance than the other methods for the three cases considered in Figure 2.

Next, we consider an example with a nonlinear RFC of the form $kx+\gamma x3$. Nonlinear RFC is considered since brain dynamics during seizures are well known to be highly nonlinear (see Stam, 2005). We vary the number of observable masses and keep the other parameters the same as in the previous example. The value of $\gamma $ is set to one for our simulations. The results of the application of different methods for this example are shown in Figure 3. Note that the performance of the SLRPM method is comparable to or better than that of the other methods for a varying number of observed variables. We also have considered a nonlinear RFC of the form $kx+\gamma x9$, and in that case too, the SLRPM method outperformed the correlation and SRPM methods (results not shown here). In Figure 4, we have plotted the performance of the methods in terms of number of samples $N$ for $p=100$. We note that the performances of the correlation, SRPM, and SLRPM methods do not improve with an increase in $N$, with the SLRPM method outperforming the correlation and SRPM methods. The correlation and SRPM methods also have comparable performance. It is worth noting that although the SLRPM method performs much better and the PM method is the worst when the sample size is small, the PM method is able to outperform all the other methods when we increase the sample size.

In the previous two examples, the networks were sparse (all masses were connected only to their immediate neighbors). For an epileptic human brain, this sparseness assumption might not remain valid (Steriade et al., 1993). Hence, we now consider relatively denser networks and evaluate the performance of the methods. In the next two examples, we consider the spring-mass network to have long-range connections. First, assume that this range is $40n$, where $n$ is a natural number (i.e., the $i$th mass); besides being connected to its immediate neighbors $i+1$ and $i-1$, it is also now connected to masses with indices $i+40$, $i-40$, $i+80$, $i-80$, and so on. We also assume a linear RFC for this example. Using Newton's second law of motion and Hooke's law, we can easily write the displacement equations of the masses for this modified spring-mass network as done for the two previous examples in equation 2.1. We can then iteratively solve these equations to generate $N$ samples of the displacement vector. We assume $N=100,000$ samples for this example. All other parameters are the same as before. Like the previous examples, the values of observable masses $p$ are taken to be 60, 100, and 150. The performance of the methods for this example is shown in Figure 5. The values of the regularization parameters are also indicated in that figure. We note that the SLRPM method is able to outperform the other methods. It is also worth noting that the performance of the SRPM method is comparable to that of the SLRPM method.

Next, assume that the range is $20n$ (i.e., the $i$th mass). Besides being connected to its immediate neighbors $i+1$ and $i-1$, it is now also connected to masses with indices $i+20$, $i-20$, $i+40$, $i-40$, and so on. We vary the number of observable masses and keep the other parameters the same as in the previous example. The performance of different methods for this example is shown in Figure 6. The SLRPM method is able to outperform the other methods. Furthermore, in Figure 7, we have plotted the performance of the methods in terms of number of samples $N$ for $p=100$. Note that the PM method is unable to estimate the observable network, resulting in large errors. Most of these errors are numerical, since the sample covariance matrix is a poor estimator of the eigenvalues of the covariance matrix. The SLRPM method performs better than the correlation and SRPM methods. It is also worth noting that the correlation and SRPM methods have comparable performance.

In the next example, we still consider relatively dense connections, but we assume the spring-mass network to have short-range connections. We assume a neighborhood of 3 (i.e., the $i$th mass), is now connected to masses with indices $i+1$, $i-1$, $i+2$, $i-2$, $i+3$, and $i-3$. We also assume $N=100,000$ samples and a linear RFC. As before, the values of observable masses $p$ are taken to be 60, 100, and 150. The performance of the methods for this example is shown in Figure 8. The values of the regularization parameters are also mentioned in that figure. We note that the SLRPM method outperforms the correlation and PM methods and has performance comparable to that of the SRPM method. We have also plotted the performance of the methods in terms of number of samples $N$ for $p=100$ in Figure 9. The PM method is not able to have any useful estimate of the observed network due to numerical errors. The SLRPM method outperforms the correlation method and has performance comparable to the SRPM method.

Although for each case, we have considered only one random subset of network to be observable, considering different random subsets of network (but fixing all the parameters) to be observable does not change the trends of the results reported here^{1} and the SLRPM method is able to outperform the other methods for different random subsets of observed variables.

The performance (in terms of error percentage) of each method for network recovery in the spring-mass model does not change under different types of noise distributions. We have tested the performance of these methods under gaussian, Poisson, and uniform distributions, and each method has similar performance under these different noise distributions.

The performance (in terms of error percentage) of each method for network recovery in the spring-mass model does not change for a change in signal-to-noise ratio up to 30 dB.

Remarks ^{2} and ^{3} are not surprising. Note that in the spring-mass model, we model the noise as an external force. So no matter how forcefully (noise variance) or in what way (noise distribution) we wobble the spring-mass system, the connectivity pattern will not change. Noise is not responsible for the connectivity pattern, and the connectivity patterns estimated by the different methods shown here are due to the methods themselves.

## 3 Epileptic Seizure Analysis from Human ECoG Recordings

### 3.1 ECoG Data Acquisition and Protocol

Continuous ECoG recordings from five patients (see Table 1 for details) with long-standing pharmaco-resistant complex partial epileptic seizures were analyzed. Recordings were performed using a standard clinical recording system (XLTEK, Natus Medical, San Carlos, CA) with a 500 Hz sampling rate. The reference channel was a strip of electrodes placed outside the dura and facing the skull at a region remote from the other grid and strip electrodes. Subdural electrode arrays were placed to confirm the hypothesized seizure focus and locate epileptogenic tissue in relation to essential cortex, thus directing surgical treatment. The decision to implant, the electrode targets, and the duration of implantation were made entirely on clinical grounds, with no input from this research study. All data acquisition was performed under protocols monitored by the Institutional Review Board of the Massachusetts General Hospital according to National Institutes of Health guidelines.

### 3.2 Preprocessing and Referencing

ECoG recordings were first low-pass-filtered at 125 Hz using a sixth-order Butterworth filter to remove high-frequency artifacts. Line frequencies 60 Hz and 120 Hz were then notch-filtered using a fourth-order Butterworth filter. Next, to reduce the signals from the reference electrode, at each time point, the average signal of all electrodes was subtracted from each electrode (Cimenser et al., 2011; Kramer et al., 2010; Nunez & Srinivasan, 2006), a process also known as common average referencing (CAR). Finally, recordings were $z$-scored (mean-variance normalization) for each channel (Varsavsky, Mareels, & Cook, 2010).

After preprocessing the ECoG data, we analyzed them using the correlation, SRPM, and SLRPM methods and used MI, CC, and EC to quantify brain connectivity. These three measures, widely used for characterizing activity in brain imaging studies, are directly applied on the estimated connectivity matrix from the methods considered. In order to calculate the measures, we used the brain connectivity toolbox (BCT; Sporns, 2010).

### 3.3 Description of Brain Connectivity Measures

#### 3.3.1 Modularity Index (MI)

#### 3.3.2 Clustering Coefficient (CC)

#### 3.3.3 Eigenvector Centrality (EC)

The largest positive eigenvalue $\kappa 1$ of $B$ is always greater than or equal to the magnitude of its most negative eigenvalue.

See Newman (2010, p. 346).$\u25a1$

There exists an eigenvector with eigenvalue $\kappa 1$ all of whose elements are nonnegative.

This follows from theorem ^{4}. See Newman (2010, p. 347), for details.$\u25a1$

Since there can be at most one eigenvector with all elements nonnegative, the eigenvector corresponding to the largest positive eigenvalue contains the relative ECs of the brain regions. The optimization problem in equation 3.5 can be solved by using the power method (Newman, 2010). The EC measure has been previously used by researchers to quantify connectivity in human brain imaging studies (Burns et al., 2014; Joyce, Laurienti, Burdette, & Hayasaka, 2010; Lohmann et al., 2010; Zuo et al., 2012).

### 3.4 Seizure Analysis

After preprocessing of the ECoG recordings, the methods are applied on 4 s nonoverlapping time windows, and brain connectivity is then quantified using the described measures. The length of the time series analyzed for each seizure was 25 minutes, with a 10 minute preictal time segment. Note that although it is difficult to characterize the preictal period of seizures, which varies from seizure to seizure within a patient and in seizures from different patients, we define the 10 minute preseizure period as our preictal period for all seizures in all patients.

The regularization parameter $\lambda $ in the SRPM method was set to 0.02, and the regularization parameters $\alpha $ and $\beta $ in the SLRPM method were set to 0.02 and 0.2, respectively, for all seizures and all patients. The SRPM and SLRPM methods were found to be robust to changes in the regularization parameters, and small changes in the values of these did not change the results and conclusions of the letter.

We first describe seizures from patient 1 in detail. Electrode locations are shown in Figure 10. As shown, there were two grids: one over the anterior temporal and fronto-parietal region and the other placed over the posterior temporal-parietal region. Three strips were placed in the frontopolar region, the subfrontal region, and the subtemporal region. Depth electrodes were inserted in the inferior frontal and the anterior and posterior temporal regions. We analyzed five clinical seizures from this patient.

The MI plots for the three methods for seizure 1 are shown in Figure 11. This seizure lasted for approximately 90 seconds (determined by the clinician). Observe that the MI estimated by the SRPM and SLRPM methods during the seizure are relatively lower than the preictal and postictal MI estimated by the same methods, respectively. Furthermore, the preictal and postictal MI in the SLRPM method are relatively higher than those in the SRPM method, respectively, and both methods have higher MI than the correlation method. Since the SLRPM method estimates brain connectivity by considering both the common observable and latent inputs, it removes considerable between-module connectivity, and the estimated brain network by this method is highly modular. Also note that the MI shown in the correlation method slightly increases during the seizure and then decreases for the first few seconds of the postictal period. These dynamics estimated by the three methods can be more clearly interpreted if we visualize the putative latent inputs to the electrodes from which brain activity was recorded. In order to have a measure for the latent inputs, we calculate the sum of the eigenvalues of the low-rank matrix estimated by the SLRPM method. This plot is shown in Figure 12. We see that the putative latent inputs are relatively low in the ictal period than the preictal period but are relatively high in the postictal period. For the first few seconds of the postictal period, these latent inputs start increasing and make the observed brain regions more correlated and consistent with a decrease in MI in the correlation method for this period. Also, since the effect of the latent inputs is relatively higher in the preictal and postictal periods than the ictal period, it becomes essential to estimate the precision matrix of the conditional statistics in order to remove the influence of these inputs on the recorded activity. SRPM and correlation methods have lower MI in preictal and postictal periods than the SLRPM method since they cannot model the latent inputs and might produce erroneous connections among the brain regions, resulting in an erroneous estimate of MI.

The low MI values in the SRPM and SLRPM methods can be explained by the CC and EC plots for the individual electrodes, which are shown in Figures 13 and 14, respectively. Note from the CC and EC plots of the SRPM and SLRPM methods that the brain regions are uniformly active in the ictal period, in contrast to the preictal and postictal periods. This relatively more uniform activity explains the lower MI shown in these two methods since during the ictal period, the entire brain network behaves as a single module. For both CC and EC measures, the SLRPM method shows relatively sparser activity than the SRPM and correlation methods during the preictal and postictal periods, which is again a consequence of the ability of the SLRPM method to model the latent inputs. The estimates by the SRPM and correlation methods in these periods might be erroneous due to the inability of these two methods to capture the effect of the latent inputs on the observed brain regions. For the first few seconds of the postictal period, the EC measure in the correlation method shows that the electrodes are relatively highly active. This phenomenon is also observed in the CC measure of the correlation method to some extent, in which we see that a relatively large number of electrodes are uniformly active. These probably are erroneous estimates due to the presence of the common inputs, as explained before. In contrast, examining the CC and EC measures in the SLRPM method during this period, we see that relatively few electrodes are active. All of these results demonstrate the superior ability of the SLRPM method to capture brain dynamics that might be the closest to the ground truth.

Plots of MI, CC, EC, and latent inputs measures from the SLRPM method for the other four seizures are shown in Figures 15, 16, 17, and 18, respectively. Each of these seizures lasted approximately 90 seconds. We note that ictal dynamics for all these seizures for each measure are very similar to each other and seizure 1 and can be explained as before.

Due to the consistency of dynamics of all five seizures in the SLRPM method, we were able to write a seizure detection algorithm for this patient. For seizure detection, we considered the latent inputs measure. Overall, we found this measure to be the most sensitive to seizures from the analysis of all seizures from all patients and hence used it for seizure detection. We used both threshold criteria, when the latent inputs become lower than a specified threshold, and duration criteria, when the duration of the latent inputs lower than the specified threshold is higher than some specified lower limit and lower than some specified higher limit. We declare that a seizure is detected when the latent inputs satisfy both of these criteria. Using these criteria, we were able to detect all five seizures with zero false positives.

Dynamics of seizures in patients 2, 3, 4, and 5 were very similar to that in patient 1 (see the appendix for detailed analysis of these patients' seizures). We also used the threshold and duration criteria on the latent inputs measure for seizure detection in these patients (the threshold and duration values for different patients were different from one another and also from those selected for patient 1) and were able to detect 5 out of 5 seizures with zero false positives in patient 2, 7 out of 7 seizures with 2 false positives in patient 3, 3 out of 3 seizures with zero false positives in patient 4, and 52 out of 54 seizures with 2 false positives in patient 5. Moreover, for patient 5, 6 seizures were not marked clinically, but showed clear signs of ictal dynamics in the latent inputs plots (and in plots of other measures not shown here) in Figure 19 (2 seizures in Figure 19a and one each in Figures 19b to 19e). This shows the potential of the SLRPM method to detect seizures that the clinician possibly missed.

### 3.5 Effect of Preprocessing

Here, we study the effect of two of the most important preprocessing steps that we applied on our ECoG recordings: CAR and $z$-scoring.

#### 3.5.1 Effect of CAR

In our CAR, we subtracted the average signal of all electrodes from each electrode at each time point. However, since the dynamics of the seizure onset electrodes might be different from those of the nononset electrodes, we tested the SLRPM method by considering only the nononset electrodes for CAR. We also tested SLRPM by considering a random subset of 40 nononset electrodes for CAR. We considered the EC measure applied on seizure 1 in patient 1 as an example to test these two CAR procedures, the results are shown in Figures 20a and 20b, respectively. For comparison, the EC plot of SLRPM method in Figure 13 is reproduced in Figure 20c in which we use all electrodes for CAR. We note that both of these CAR approaches have minimal effect on the preictal, ictal, and postictal dynamics, and the dynamics look similar to that in Figure 20c.

#### 3.5.2 Effect of normalization

For our analysis, we $z$-scored each channel for the entire 25 minute time segment to have a relatively less biased estimation of mean and standard deviation. However, some researchers also use $z$-scoring in each time window rather than the entire data (Varsavsky et al., 2010). To test the effect of normalization in each time window on SLRPM, we considered the EC measure applied on seizure 1 in patient 1 as an example. Figure 21a shows the result of this normalization procedure. For comparison, the EC plot of SLRPM method in Figure 13 is reproduced in Figure 21b in which normalization was applied for the entire time segment. We note that the preictal dynamics (also see Figure 22 for the plot of only the preictal period of Figure 21a) in Figure 21a is significantly different from that in Figure 21b and the ictal and postictal dynamics are relatively more similar for these two normalization procedures. Moreover, we have applied $z$-scoring only on the 10 minute preictal recording; the result is shown in Figure 21c. Note that although the values of the EC measures for the electrodes are different for the preictal period in Figures 21b and 21c, electrodes that were relatively more active than others can be seen in both. However, both of these are markedly different from the EC measures for the electrodes in the preictal period in Figure 21a. This shows that two different normalization methods can produce different results, and care needs to be taken in interpreting the results.

## 4 Discussion

Analyzing the results from the measures applied on the SLRPM method, we see that ictal dynamics is relatively more consistent than the preictal dynamics across seizures within a patient. Researchers (Kramer et al., 2010; Ponten et al., 2007; Schindler et al., 2008; Vega-Zelaya et al., 2015) have previously reported an increase in ACC after seizure onset. However, in our analysis, we do not consider ACC since we found that in the preictal period, some electrodes are relatively more active than others, and averaging CC across electrodes in this period and comparing the result to that in the ictal period obscures the dynamics for individual electrodes. We also observe a relative decrease in MI after seizure onset in patients 1 to 4, and this result is consistent with findings from Vega-Zelaya et al. (2015). Similar to the authors in Burns et al. (2014), we also observed change in brain connectivity from the EC measure. However, we were not able to locate the seizure onset zones from the EC measure applied on the SLRPM method. The SLRPM method estimates the conditional statistics of the observed brain regions given the latent brain regions hence removing the effect of the latent inputs on the connectivity of the measured brain regions. We found the estimated latent inputs by the SLRPM method to be the most sensitive measure to the seizures in the ECoG recordings from the five patients we analyzed, and using this measure, we were able to detect seizures with very high accuracy (72 out of 74 seizures) and very low false positives (4). The latent inputs measure performed better than both correlation and SRPM methods for seizure detection (for correlation and SRPM methods, we used the MI measure for seizure detection). A comprehensive comparison of our seizure detection algorithm with other seizure detection methods Jerger et al., 2001; e.g., frequency-based methods) will be explored in a future paper to place our seizure detection method in the context of other algorithms. There is also evidence of decreased cortical-subcortical interaction during epileptic seizures (Blumenfeld et al., 2012), which explains the relatively decreased latent inputs during the ictal period estimated by the SLRPM method. There was a noticeable offset between the clinical onset times (green lines) and the onset of the uniform region of CC in Figure 25 (best seen in panel c). A similar offset in MI is apparent in Figure 24c and in the putative latent input in Figure 27c. This suggests that these measures may be useful for automating the detections of seizure onsets. Furthermore, we also observed variability in preictal activity across seizures in each patient from our measures. Quantification of this heterogeneity and detailed interpretation of preseizure brain activity for patients are published elsewhere (Das, 2018; Das, Cash, & Sejnowski, 2018).

An alternative way of estimating the internal states or latent inputs of brain activity from a few ECoG electrode recordings is using the observability framework from control systems theory (Liu, Slotine, & Barabasi, 2013; Meyer-Base et al., 2017; Whalen, Brennan, Sauer, & Schiff, 2015). If the complete internal state dynamics of the brain can be reconstructed from the available ECoG recordings, then we will have an observable brain (Aguirre, Portes, & Letellier, 2017). Controlling the observable brain regions responsible for seizure initiation will be important to gain significant insights into its evolution, thus enabling us to find better treatments of epilepsy. However, this problem is particularly confounded by the highly nonlinear, nonstationary, and complex nature of the human brain, for which the ground-truth topology is also unknown. Moreover, the nonlinear, long-range connections considered in the spring-mass network in this letter are already too complicated to carry out a theoretical analysis from a controllability and observability point of view. Nevertheless, we have carried out an empirical analysis of such networks as considered previously (Huth et al., 2016), and this has helped us validate the performance of the methods considered in this letter and the effect of internal states or latent inputs on the inferred connectivity by such methods. However, epileptic seizure generation and propagation are still poorly understood, and clinically relevant signal processing methods are scarce. Combining graph theory–based methods with the theory of observability and controllability of human brain networks seems to be the future of tackling such a highly complex problem (Whalen et al., 2015). This will be crucial in revealing the underlying neural circuits in epileptic seizures and more broadly will have significant impact on our understanding of the human brain in pathology.

## Appendix

### A.1 Medication Information of Patients

Medication information for patient 1 is shown in Table 2. Seizure 1 corresponds to day 7, seizures 2 and 3 correspond to day 8, and seizures 4 and 5 correspond to day 9. Medication information for patient 2 is shown in Table 3. Seizures 1 and 2 correspond to day 2, seizure 3 corresponds to day 3, and seizures 4 and 5 correspond to day 8. Medication information for patient 3 is shown in Table 4. Seizure 1 corresponds to day 2; seizures 2, 3, 4, and 5 correspond to day 7; and seizures 6 and 7 correspond to day 8. Medication information for patient 4 was not available. Patient 5's hospital medications included phenobarbital, carbamazepine, ativan, and vimpat.

Day . | Medication . |
---|---|

On admission | Trileptal 300 tid + 150 mg OD, clonazepam 0.25 mg tds |

Day 1 | Trileptal 300 bid; clonazepam 0.25 mg tds |

Day 2 | Trileptal 300 od; clonazepam 0.25 mg tds |

Day 3 | Clonazepam 0.25 mg tds |

Day 4 | Clonazepam 0.25 mg bid |

Day 5 | Clonazepam 0.25 mg od |

Day $6--9$ | No AEDs |

Day $10--11$ and on discharge | Trileptal 300 tid + 150 mg OD, clonazepam 0.25 mg tds |

Day . | Medication . |
---|---|

On admission | Trileptal 300 tid + 150 mg OD, clonazepam 0.25 mg tds |

Day 1 | Trileptal 300 bid; clonazepam 0.25 mg tds |

Day 2 | Trileptal 300 od; clonazepam 0.25 mg tds |

Day 3 | Clonazepam 0.25 mg tds |

Day 4 | Clonazepam 0.25 mg bid |

Day 5 | Clonazepam 0.25 mg od |

Day $6--9$ | No AEDs |

Day $10--11$ and on discharge | Trileptal 300 tid + 150 mg OD, clonazepam 0.25 mg tds |

Day . | Medication . |
---|---|

Home doses | Levetiracetam $1500--1500$, oxcarbazepine $600--600$, clonazepam $0.25--0.25$ |

Day 1 | Levetiracetam $1500-1500$, oxcarbazepine $600--600$, clonazepam $0.25--0.25$ |

Day 2 | Levetiracetam $1500--1500$, oxcarbazepine $300--300$, clonazepam $0.25--0.25$ |

Day 3 | Levetiracetam $1500--1500$, oxcarbazepine $300--600$, clonazepam $0.25--0.25$ |

Day 4 | Levetiracetam $1500--1500$, oxcarbazepine $600--600$, clonazepam $0.25--0.25$ |

Day 5 | Levetiracetam $1500--750$, oxcarbazepine $600--300$, clonazepam $0.25--0.25$ |

Day 6 | Levetiracetam $750--750$, oxcarbazepine $300--0$, clonazepam stopped |

Day 7 | Levetiracetam stopped, oxcarbazepine stopped, clonazepam stopped |

Day 8 | Levetiracetam $750--750$, oxcarbazepine stopped, clonazepam stopped |

Day 9 | Levetiracetam $1500--1500$, oxcarbazepine $600--600$, clonazepam $0.25--0.25$ |

Day . | Medication . |
---|---|

Home doses | Levetiracetam $1500--1500$, oxcarbazepine $600--600$, clonazepam $0.25--0.25$ |

Day 1 | Levetiracetam $1500-1500$, oxcarbazepine $600--600$, clonazepam $0.25--0.25$ |

Day 2 | Levetiracetam $1500--1500$, oxcarbazepine $300--300$, clonazepam $0.25--0.25$ |

Day 3 | Levetiracetam $1500--1500$, oxcarbazepine $300--600$, clonazepam $0.25--0.25$ |

Day 4 | Levetiracetam $1500--1500$, oxcarbazepine $600--600$, clonazepam $0.25--0.25$ |

Day 5 | Levetiracetam $1500--750$, oxcarbazepine $600--300$, clonazepam $0.25--0.25$ |

Day 6 | Levetiracetam $750--750$, oxcarbazepine $300--0$, clonazepam stopped |

Day 7 | Levetiracetam stopped, oxcarbazepine stopped, clonazepam stopped |

Day 8 | Levetiracetam $750--750$, oxcarbazepine stopped, clonazepam stopped |

Day 9 | Levetiracetam $1500--1500$, oxcarbazepine $600--600$, clonazepam $0.25--0.25$ |

Day . | Medication . |
---|---|

Home doses | Keppra $1500--1500$, lamotrigine $200--200$, Zonegran $300--300$ |

Day 1 | Keppra $1500--1500$, lamotrigine $200--200$, Zonegran $300--300$ |

Day 2 | Keppra $750--750$, lamotrigine $100--100$, Zonegran $200--200$ |

Day 3 | Keppra $375--375$, lamotrigine $50--50$, Zonegran $100--100$ |

Day $4--8$ | Off medications |

Day . | Medication . |
---|---|

Home doses | Keppra $1500--1500$, lamotrigine $200--200$, Zonegran $300--300$ |

Day 1 | Keppra $1500--1500$, lamotrigine $200--200$, Zonegran $300--300$ |

Day 2 | Keppra $750--750$, lamotrigine $100--100$, Zonegran $200--200$ |

Day 3 | Keppra $375--375$, lamotrigine $50--50$, Zonegran $100--100$ |

Day $4--8$ | Off medications |

### A.2 Seizure Characteristics in Patient 2

Electrode locations are shown in Figure 23. As shown, depth electrodes were inserted in the right and left anterior temporal regions, right and left posterior temporal regions, right and left cingulate regions, and right and left subfrontal regions. We analyzed five clinical seizures from this patient.

Plots of MI, CC, EC, and latent inputs measures from the SLRPM method for five seizures of this patient are shown in Figures 24, 25, 26, and 27, respectively. From the MI plots of these seizures, ictal dynamics is visible, to varying extent, by the relative lowering of MI for the first few seconds of seizure onset, and this phenomenon is similar to the ictal dynamics seen in patient 1. A relatively large number of electrodes are uniformly active for the first few seconds of seizure onset in each of the CC and EC plots, which is similar to those observed in patient 1. It is also noteworthy that the number of preictal active electrodes is relatively sparser in the EC plots than those in the corresponding CC plots. We also observe the relative lessening of the putative latent inputs, to varying extent, after seizure onset, a phenomenon also observed in patient 1.

### A.3 Seizure Characteristics in Patient 3

Electrode locations are shown in Figure 28. As shown, depth electrodes were inserted in the right and left anterior temporal regions, right and left posterior temporal regions, right and left cingulate regions, right and left subfrontal regions, and right and left posterior frontal regions. We analyzed seven clinical seizures from this patient.

### A.4 Seizure Characteristics in Patient 4

Electrode locations are shown in Figure 33. As shown, depth electrodes were inserted in the anterior temporal, posterior temporal, subfrontal, cingulate, and parietal regions. We analyzed three clinical seizures from this patient.

Plots of MI, CC, EC, and latent inputs measures from the SLRPM method for three seizures of this patient are shown in Figures 34, 35, 36, and 37, respectively. Ictal dynamics is clearly visible from these plots and similar to that in patients 1, 2, and 3. Large variation in MI and latent inputs plots is also observed starting from the 700 s mark and lasting for a few seconds in seizure 3 and seizure 2 to some extent. From all plots for seizure 1, we also see that the duration of ictal dynamics is relatively large as compared to those for the other two seizures.

### A.5 Seizure Characteristics in Patient 5

Electrode locations are shown in Figure 38. As shown, there were two grids, one placed over the posterior inferior frontal, posterior inferior parietal, superior temporal, and superior posterior temporal regions and the other covered the posterior inferior frontal, inferior parietal, and superior parietal lobe regions. Two strips were placed over the subtemporal and parietal regions. There were also depth electrodes placed in the anterior temporal lobe, the posterior temporal lobe, the superior parietal lobe, and encephalomalacia noted on MRI, and through the posterior frontal lobe to reach the cingulate gyrus. We analyzed 54 clinical seizures from this patient.

Plots of MI, CC, EC, and latent inputs measures from the SLRPM method for four seizures of this patient are shown in Figures 39, 40, 41, and 42, respectively. We note that ictal dynamics is not obvious from the MI plots, but can be seen, to a varying extent, in the CC plots where a relatively large number of electrodes become uniformly active. Clear change in preictal dynamics in the EC plots also demonstrates the capability of the SLRPM method to capture ictal dynamics. The latent inputs measure appears to be the most sensitive among all measures in this patient. Ictal dynamics is very clear by the relative lowering of the putative latent inputs followed by a sharp increase.

## Note

^{1}

Figures 2, 3, 5, 6, and 8 also appear in Das (2018). The figures in this paper and the figures in Das (2018) were generated in two different versions of Matlab, and they are slightly different from each other due to numerical precision (floating-point arithmetic) issues. But the trends of the results are the same in both cases.

## Acknowledgments

This research was supported by the Howard Hughes Medical Institute, the Swartz Foundation, the NIH under grant R01EB009282, and the Office of Naval Research under grant N000141310672. We also thank David J. C. MacKay for interesting discussions on the spring-mass model, Mark A. Kramer for discussions on the preprocessing steps in the ECoG recordings, Giovanni Piantoni for his help in using the FieldTrip toolbox, and Aaron L. Sampson for discussions on epilepsy.