## 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.

However, a major drawback of the PM method is that under the regime of a relatively small number of time samples, the sample covariance matrix might not be directly invertible. Even if there are sufficient numbers of samples such that the sample covariance matrix is invertible, the estimated PM might produce large errors in connections in a given brain network since the sample covariance matrix is a poor estimator of the eigenvalues of the covariance matrix (Das et al., 2017; see the simulation results in section 2). Hence, researchers have proposed using a sparse-regularized precision matrix (SRPM) to calculate the pairwise partial correlations (Friedman, Hastie, & Tibshirani, 2008). The SRPM can be estimated by solving the following $L1$ regularized optimization problem for $X$,
$argminXs.t.X≻0-logdet(X)+tr(SX)+λ∥X∥1,$
(1.1)
where $λ$ is the regularization parameter balancing the error in the maximum likelihood estimate (MLE) of the precision matrix and the sparsity (the MLE of the precision matrix is the inverse of the sample covariance matrix according to the invariance principle) and $S$ is the sample covariance matrix. The effect of the $L1$ regularization term $λ∥X∥1$ is twofold. First, it helps make the estimation of the precision matrix problem well posed since direct inversion of the sample covariance matrix is ill posed when finite numbers of time samples are available. Second, since it promotes sparsity, it can incorporate this assumption while estimating the human brain network, which is known to be sparse (Laughlin & Sejnowski, 2003; Steriade, McCormick, & Sejnowski, 1993). Observe that the optimization problem in equation 1.1 is a convex optimization problem, and several algorithms (Banerjee, El Ghaoui, & d'Aspremont, 2008; Friedman et al., 2008; Hsieh, Sustik, Dhillon, & Ravikumar, 2011; Hsieh, Sustik, Dhillon, Ravikumar, & Poldrack, 2013; Oztoprak, Nocedal, Rennie, & Olsen, 2012; Rothman, Bickel, Levina, & Zhu, 2008; Scheinberg, Ma, & Goldfarb, 2010; Yuan & Lin, 2007) have been proposed for solving this optimization problem. For the analysis in our letter, we use the QUIC algorithm (Hsieh et al., 2011, 2013) to estimate the SRPM for its relatively fast computation time.

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

Chandrasekaran, Parrilo, and Willsky (2012) recently proposed using a sparse-plus-latent-regularized precision matrix (SLRPM) when there are unobserved or latent regions interacting with the observed regions, which is typical of brain imaging. 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. Here we assume that the observed and latent variables jointly follow a multivariate gaussian distribution. We also assume that we have $p$ number of observed variables and $h$ number of latent variables. Let the covariance matrix $Σ$ of this joint distribution be represented by
$Σ=ΣPΣP,HΣH,PΣH,$
(1.2)
where $ΣP$ is of dimension $p×p$ representing the covariances of the observed variables, $ΣP,H$ is of dimension $p×h$ representing the covariances between the observed and latent variables, and $ΣH$ is of dimension $h×h$ representing the covariances of the latent variables. Let the inverse $Σ-1$ of this covariance matrix be represented by
$Σ-1≜K=KPKP,HKH,PKH,$
(1.3)
where $KP$ is of dimension $p×p$ representing the dependencies of the observed variables, $KP,H$ is of dimension $p×h$ representing the dependencies between the observed and latent variables, and $KH$ is of dimension $h×h$ representing the dependencies of the latent variables. The marginal precision matrix $ΣP-1$ of the observed variables can be given by the Schur complement with regard to $KH$ as in Chandrasekaran et al. (2012):
$K˜P≜ΣP-1=KP-KP,HKH-1KH,P.$
(1.4)
The matrix $KP$ is the precision matrix of the conditional statistics of the observed variables given the latent variables. Hence, if the underlying brain network is sparse, then $KP$ is a sparse matrix. The matrix $KP,HKH-1KH,P$ summarizes the effect of marginalization over the latent variables. Note that the rank of $KP,HKH-1KH,P$ can be at most $h$, and, hence, it will be of small rank if the number of latent variables is relatively small. Thus, the marginal precision matrix $K˜P$ is not sparse in general due to the presence of the term $KP,HKH-1KH,P$, and the SRPM method for the observed variables might not be able to correctly identify the brain connectivity. Since from equation 1.4, the marginal precision matrix is decomposed as a sum of a sparse matrix and a low-rank matrix, we impose the $L1$ norm on the precision matrix of the conditional statistics (also known as the SLRPM) and nuclear norm on the low-rank matrix. By solving the following regularized optimization problem, $X$ provides an estimate of the sparse matrix $KP$, and $L$ provides an estimate of the low-rank matrix $KP,HKH-1KH,P$,
$argminX,Ls.tX-L≻0,L⪰0-logdet(X-L)+tr(S(X-L))+α∥X∥1+βtr(L),$
(1.5)
where $α$ and $β$ are the regularization parameters balancing the error in the likelihood and the sparse and low-rank terms and $S$ is the sample covariance matrix. The $L1$ regularization term $α∥X∥1$ imposes sparsity on the underlying brain connectivity, and the trace or nuclear norm regularization term $βtr(L)$ imposes low rankness on the common inputs from the latent or unobserved brain regions. Furthermore, these regularizations make the optimization problem well behaved when we have a finite number of samples. The optimization problem in equation 1.5 is a convex problem, and we use the alternating direction method of multipliers (ADMM; Ma, Xue, & Zou, 2013) to estimate the SLRPM for our analysis. Note that the SRPM and SLRPM have to be normalized to obtain the partial correlations (Whittaker, 1990).

### 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

Consider the spring-mass model shown in Figure 1. We assume that we can measure the displacements of the masses in the model. Also assume there are 200 equal masses $m$ connected in cascade via springs, that is, the $i$th mass is connected to its neighbors $i-1$ and $i+1$ via springs. First, consider each spring to have a restoring force characteristic (RFC) of the form $kx$, where $k$ is the spring constant and $x$ is the amount by which the spring is displaced from its relaxed position. Note that the RFC is linear in this case. The leftmost and the rightmost springs are connected to a rigid wall. For our simulations (hereafter, all units are in mks), we use $m=0.1$ and $k=1$. We assume that this spring-mass model is subject to thermal perturbation. In addition, we assume that each mass is subject to external force, and we model this external force as a white noise process with variance $σ2$. For our simulation, we use $σ2=0.000025$. We denote the displacements of the masses as $x1,x2,…,x200$ and the external forces associated with each of them as $w1,w2,…,w200$, respectively. Using Newton's second law of motion and Hooke's law, we can write the displacement equations of the masses as
$mx¨1=-kx1+k(x2-x1)+w1mx¨2=-k(x2-x1)+k(x3-x2)+w2⋮mx¨200=-k(x200-x199)-kx200+w200.$
(2.1)
In matrix-vector form, the above set of equations can be written as
$mx¨1x¨2x¨3x¨4⋮x¨200=k-2100⋯01-210⋯001-21⋯0001-2⋯0⋮⋮⋮⋮⋱⋮000⋯1-2x1x2x3x4⋮x200+w1w2w3w4⋮w200,$
(2.2)
where $x1,x2,…,x200$ measure the displacements relative to the equilibrium positioning.
Figure 1:

The spring-mass model. The black dots denote continuation of springs and masses.

Figure 1:

The spring-mass model. The black dots denote continuation of springs and masses.

In compact form, equation 2.2 can be represented as
$mx¨=kCx+w,$
(2.3)
where $w$ is the white noise vector and $C$ is the connectivity (or ground-truth) matrix of the entire spring-mass network. Using a second-order approximation of the double derivative on the left-hand side in equation 2.3, we have
$mx(t+h)-2x(t)+x(t-h)h2=kCx(t)+w(t),$
(2.4)
where $h$ is the step size and $t$ denotes the time instant. We use a random initialization of $x$ with variance $σx2=0.000001$. We then solve equation 2.4 repeatedly to generate $N$ samples of the displacement vector $x$ corresponding to $N$ consecutive time points. For our simulation, we choose $h=0.0007$ and $N=50,000$. We can use these $N$ samples to form the sample covariance matrix of the entire spring-mass network.

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.

Figure 2:

Results of the estimation methods for the spring-mass model of 200 masses when only part of the network is observable and a linear RFC is considered. Connections in the observed network (ground-truth matrix) are shown in black dots. For the methods, correct identification of a connection is shown with black dots and incorrect identification of a connection with red dots. (a) $p=60$. (b) $p=100$. (c) $p=150$.

Figure 2:

Results of the estimation methods for the spring-mass model of 200 masses when only part of the network is observable and a linear RFC is considered. Connections in the observed network (ground-truth matrix) are shown in black dots. For the methods, correct identification of a connection is shown with black dots and incorrect identification of a connection with red dots. (a) $p=60$. (b) $p=100$. (c) $p=150$.

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 $λ$ for the SRPM method and the values of the regularization parameters $α$ and $β$ 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+γ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 $γ$ 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+γ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.

Figure 3:

Results of the estimation methods for the spring-mass model of 200 masses when only part of the network is observable and a nonlinear RFC is considered. Connections in the observed network (ground-truth matrix) are shown in black dots. For the methods, correct identification of a connection is shown in black dots and incorrect identification of a connection is shown in red dots. (a) $p=60$. (b) $p=100$. (c) $p=150$.

Figure 3:

Results of the estimation methods for the spring-mass model of 200 masses when only part of the network is observable and a nonlinear RFC is considered. Connections in the observed network (ground-truth matrix) are shown in black dots. For the methods, correct identification of a connection is shown in black dots and incorrect identification of a connection is shown in red dots. (a) $p=60$. (b) $p=100$. (c) $p=150$.

Figure 4:

Plot of error percentage of the estimation methods versus $N$ with $p=100$ and when a nonlinear RFC is considered.

Figure 4:

Plot of error percentage of the estimation methods versus $N$ with $p=100$ and when a nonlinear RFC is considered.

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.

Figure 5:

Results of the estimation methods for the spring-mass model of 200 masses when only part of the network is observable and there are long-range connections with range $40n$. Connections in the observed network (ground-truth matrix) are shown in black dots. For the methods, correct identification of a connection is shown in black dots and incorrect identification of a connection in red dots. (a) $p=60$. (b) $p=100$. (c) $p=150$.

Figure 5:

Results of the estimation methods for the spring-mass model of 200 masses when only part of the network is observable and there are long-range connections with range $40n$. Connections in the observed network (ground-truth matrix) are shown in black dots. For the methods, correct identification of a connection is shown in black dots and incorrect identification of a connection in red dots. (a) $p=60$. (b) $p=100$. (c) $p=150$.

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.

Figure 6:

Results of the estimation methods for the spring-mass model of 200 masses when only part of the network is observable and there are long-range connections with range $20n$. Connections in the observed network (ground-truth matrix) are shown in black dots. For the methods, correct identification of a connection is shown in black dots and incorrect identification of a connection in red dots. (a) $p=60$. (b) $p=100$. (c) $p=150$.

Figure 6:

Results of the estimation methods for the spring-mass model of 200 masses when only part of the network is observable and there are long-range connections with range $20n$. Connections in the observed network (ground-truth matrix) are shown in black dots. For the methods, correct identification of a connection is shown in black dots and incorrect identification of a connection in red dots. (a) $p=60$. (b) $p=100$. (c) $p=150$.

Figure 7:

Plot of error percentages of the estimation methods versus $N$ with $p=100$ and a densely connected network with long-range connections of range $20n$.

Figure 7:

Plot of error percentages of the estimation methods versus $N$ with $p=100$ and a densely connected network with long-range connections of range $20n$.

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.

Figure 8:

Results of the estimation methods for the spring-mass model of 200 masses when only part of the network is observable and there are short-range connections with a neighborhood of 3. Connections in the observed network (ground-truth matrix) are shown in black dots. For the methods, correct identification of a connection is shown in black dots and incorrect identification in red dots. (a) $p=60$. (b) $p=100$. (c) $p=150$.

Figure 8:

Results of the estimation methods for the spring-mass model of 200 masses when only part of the network is observable and there are short-range connections with a neighborhood of 3. Connections in the observed network (ground-truth matrix) are shown in black dots. For the methods, correct identification of a connection is shown in black dots and incorrect identification in red dots. (a) $p=60$. (b) $p=100$. (c) $p=150$.

Figure 9:

Plot of error percentage of the estimation methods versus $N$ with $p=100$ and densely connected network with short-range connections of neighborhood 3.

Figure 9:

Plot of error percentage of the estimation methods versus $N$ with $p=100$ and densely connected network with short-range connections of neighborhood 3.

Remark 1.

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 here1 and the SLRPM method is able to outperform the other methods for different random subsets of observed variables.

Remark 2.

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.

Remark 3.

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.

Table 1:
Patient Details.
PatientAgeSex (Male (M)/Female (F))
46
55
29
45
25
PatientAgeSex (Male (M)/Female (F))
46
55
29
45
25

### 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)

For clustering of brain regions, we use MI (Blondel et al., 2008; Brandes et al., 2008; Newman & Girvan, 2004; Newman, 2006; Reichardt & Bornholdt, 2006; Ronhovde & Nussinov, 2009; Sporns, 2010; Sun et al., 2009). MI is defined as
$Q=12M∑i,jBij-sisj2Mδ(σi,σj),$
(3.1)
where $Bij$ denotes the strength of connectivity (obtained from the methods) between brain regions $i$ and $j$, $si=∑jBij$ denotes the sum of connectivity strengths between brain region $i$ and the rest of the brain regions, $σi$ denotes the cluster to which brain region $i$ belongs $M=12∑ijBij$, and $δ(σi,σj)$ is 1 if $σi=σj$ and 0 otherwise. Note that MI measures the fraction of the strengths of the connectivities belonging to the same cluster minus the probability of the strengths of random connectivities belonging to the same cluster. Hence, in order to have positive MI, the strengths of the connectivities belonging to the same cluster have to be better than random. Extraction of brain networks is based on MI optimization, and we use the Louvain method for community detection (LMCD) (Blondel et al., 2008; Brandes et al., 2008; Newman & Girvan, 2004; Newman, 2006; Reichardt & Bornholdt, 2006; Ronhovde & Nussinov, 2009; Sporns, 2010; Sun et al., 2009), which is a fully automatic method, for this optimization. In other words, LMCD maximizes the strengths of connectivities within clusters and minimizes the strengths of connectivities between clusters. Clustering of brain regions based on MI optimization is a widely used method in human brain imaging data analysis (Bassett et al., 2010; 2011; Bruno et al., 2017; Cole, Bassett, Power, Braver, & Petersen, 2014; Meunier, Lambiotte, & Bullmore, 2010; Rubinov & Sporns, 2010; 2011; Sporns, 2011; Zuo et al., 2012).

#### 3.3.2  Clustering Coefficient (CC)

The CC (Onnela, Saramaki, Kertesz, & Kaski, 2005; Saramaki, Kivela, Onnela, Kaski, & Kertesz, 2007) of a brain region in a network quantifies how well its neighboring brain regions are connected. CC $Ci$ of the $i$th brain region is defined as
$Ci=1ki(ki-1)∑j,h(BijBjhBhi)1/3$
(3.2)
where $ki$ is the number of neighboring brain regions of the $i$th brain region and $Bij$ is the strength of connectivity (obtained from the methods) between brain regions $i$ and $j$. The average CC (ACC) $C$ of an entire brain network can be defined as
$C=1M∑iCi,$
(3.3)
where $M$ is the total number of observed brain regions. CC is also a widely used measure in human brain imaging data analysis (Bruno et al., 2017; Kramer et al., 2010; Ponten et al., 2007; Rubinov & Sporns, 2010; Schindler et al., 2008; Vega-Zelaya et al., 2015; Wang, Ghumare, Vandenberghe, & Dupont, 2017).

#### 3.3.3  Eigenvector Centrality (EC)

The EC (Newman, 2010) is a measure of the influence or importance of a brain region in the entire brain network. This is based on the concept that a brain region more strongly connected to high influential brain regions will have relatively higher EC. Mathematically, the relative EC $ei$ of the $i$th brain region can be written as
$ei=1κ∑jBijej,$
(3.4)
where $Bij$ is the strength of connectivity (obtained from the methods) between brain regions $i$ and $j$ and $ej$ is the relative EC of the $j$th brain region. In matrix-vector notation, the above set of equations can be compactly written as an eigenvector equation,
$Be=κe,$
(3.5)
where $B$ is the estimated connectivity matrix from the methods and $e$ is its eigenvector, which contains the relative influences of the brain regions. Hence, the relative ECs can be found by solving the eigenvector equation in equation 3.5. But note that equation 3.5 has multiple solutions. The relative ECs are always assumed to be nonnegative. The following two theorems help us find EC from equation 3.5. These theorems are valid for an estimated connectivity matrix $B$ that is real, symmetric, and with all elements nonnegative (we take the absolute values of the elements of the estimated connectivity matrix to make them nonnegative).
Theorem 1.

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

Proof.

See Newman (2010, p. 346).$□$

Theorem 2.

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

Proof.

This follows from theorem 4. See Newman (2010, p. 347), for details.$□$

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 $λ$ in the SRPM method was set to 0.02, and the regularization parameters $α$ and $β$ 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.

Figure 10:

Electrode locations in patient 1.

Figure 10:

Electrode locations in patient 1.

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.

Figure 11:

MI plots for the three methods for seizure 1 in patient 1. The green line denotes the seizure onset time.

Figure 11:

MI plots for the three methods for seizure 1 in patient 1. The green line denotes the seizure onset time.

Figure 12:

Plot of the latent inputs estimated by the SLRPM method for seizure 1 in patient 1. The green line denotes the seizure onset time.

Figure 12:

Plot of the latent inputs estimated by the SLRPM method for seizure 1 in patient 1. The green line denotes the seizure onset time.

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.

Figure 13:

CC plots for the three methods for seizure 1 in patient 1. The green line denotes the seizure onset time.

Figure 13:

CC plots for the three methods for seizure 1 in patient 1. The green line denotes the seizure onset time.

Figure 14:

EC plots for the three methods for seizure 1 in patient 1. The green line denotes the seizure onset time.

Figure 14:

EC plots for the three methods for seizure 1 in patient 1. The green line denotes the seizure onset time.

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.

Figure 15:

MI plots for the SLRPM method for seizures 2, 3, 4, and 5 (shown in panels a, b, c, and d, respectively) in patient 1. The green lines in all plots denote the seizure onset time.

Figure 15:

MI plots for the SLRPM method for seizures 2, 3, 4, and 5 (shown in panels a, b, c, and d, respectively) in patient 1. The green lines in all plots denote the seizure onset time.

Figure 16:

CC plots for the SLRPM method for seizures 2, 3, 4, and 5 (shown in panels a, b, c, and d, respectively) in patient 1. The green lines in all plots denote the seizure onset time.

Figure 16:

CC plots for the SLRPM method for seizures 2, 3, 4, and 5 (shown in panels a, b, c, and d, respectively) in patient 1. The green lines in all plots denote the seizure onset time.

Figure 17:

EC plots for the SLRPM method for seizures 2, 3, 4, and 5 (shown in panels a, b, c, and d, respectively) in patient 1. The green lines in all plots denote the seizure onset time.

Figure 17:

EC plots for the SLRPM method for seizures 2, 3, 4, and 5 (shown in panels a, b, c, and d, respectively) in patient 1. The green lines in all plots denote the seizure onset time.

Figure 18:

Plots of the latent inputs estimated by the SLRPM method for seizures 2, 3, 4, and 5 (shown in panels a, b, c, and d, respectively) in patient 1. The green lines in all plots denote the seizure onset time.

Figure 18:

Plots of the latent inputs estimated by the SLRPM method for seizures 2, 3, 4, and 5 (shown in panels a, b, c, and d, respectively) in patient 1. The green lines in all plots denote the seizure onset time.

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.

Figure 19:

Plots of the latent inputs estimated by the SLRPM method for five seizures in patient 5 in which ictal dynamics is also observed at time instants (two in panel a and one each in panels b to e), which were not classified by the clinician as seizure onset times. The green lines in all plots denote the seizure onset time.

Figure 19:

Plots of the latent inputs estimated by the SLRPM method for five seizures in patient 5 in which ictal dynamics is also observed at time instants (two in panel a and one each in panels b to e), which were not classified by the clinician as seizure onset times. The green lines in all plots denote the seizure onset time.

### 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.

Figure 20:

Analysis of effect of CAR on SLRPM method for seizure 1 in patient 1. The green lines in all plots denote the seizure onset time. (a) Nononset electrodes are used for CAR. (b) A random subset of 40 nononset electrodes is used for CAR. (c) All electrodes are used for CAR.

Figure 20:

Analysis of effect of CAR on SLRPM method for seizure 1 in patient 1. The green lines in all plots denote the seizure onset time. (a) Nononset electrodes are used for CAR. (b) A random subset of 40 nononset electrodes is used for CAR. (c) All electrodes are used for CAR.

#### 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.

Figure 21:

Analysis of the effect of normalization on SLRPM method for seizure 1 in patient 1. The green lines in panels a and b denote the seizure onset time. (a) Window-wise $z$-scoring. (b) $Z$-scoring the entire time segment. (c) $Z$-scoring only the preictal time segment.

Figure 21:

Analysis of the effect of normalization on SLRPM method for seizure 1 in patient 1. The green lines in panels a and b denote the seizure onset time. (a) Window-wise $z$-scoring. (b) $Z$-scoring the entire time segment. (c) $Z$-scoring only the preictal time segment.

Figure 22:

Shown here is the preictal period of Figure 21a for better visualization of the active electrodes.

Figure 22:

Shown here is the preictal period of Figure 21a for better visualization of the active electrodes.

## 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.

Table 2:
Medication Information for Patient 1.
DayMedication
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
DayMedication
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
Table 3:
Medication Information for Patient 2.
DayMedication
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$
DayMedication
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$
Table 4:
Medication Information for Patient 3.
DayMedication
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
DayMedication
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.

Figure 23:

Electrode locations in patient 2. (a) Left hemisphere. (b) Right hemisphere.

Figure 23:

Electrode locations in patient 2. (a) Left hemisphere. (b) Right hemisphere.

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.

Figure 24:

MI plots for the SLRPM method for seizures 1, 2, 3, 4, and 5, shown in panels a, b, c, d, and e, respectively, in patient 2. Green lines in all plots denote the seizure onset time.

Figure 24:

MI plots for the SLRPM method for seizures 1, 2, 3, 4, and 5, shown in panels a, b, c, d, and e, respectively, in patient 2. Green lines in all plots denote the seizure onset time.

Figure 25:

CC plots for the SLRPM method for seizures 1, 2, 3, 4, and 5, shown in panels a, b, c, d and e, respectively, in patient 2. Green lines in all plots denote the seizure onset time.

Figure 25:

CC plots for the SLRPM method for seizures 1, 2, 3, 4, and 5, shown in panels a, b, c, d and e, respectively, in patient 2. Green lines in all plots denote the seizure onset time.

Figure 26:

EC plots for the SLRPM method for seizures 1, 2, 3, 4, and 5, shown in panels a, b, c, d, and e, respectively, in patient 2. Green lines in all plots denote the seizure onset time.

Figure 26:

EC plots for the SLRPM method for seizures 1, 2, 3, 4, and 5, shown in panels a, b, c, d, and e, respectively, in patient 2. Green lines in all plots denote the seizure onset time.

Figure 27:

Plots of the latent inputs estimated by the SLRPM method for seizures 1, 2, 3, 4, and 5, shown in panels a, b, c, d, and e, respectively, in patient 2. Green lines in all plots denote the seizure onset time.

Figure 27:

Plots of the latent inputs estimated by the SLRPM method for seizures 1, 2, 3, 4, and 5, shown in panels a, b, c, d, and e, respectively, in patient 2. Green lines in all plots denote the seizure onset time.

### 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.

Figure 28:

Electrode locations in patient 3. (a) Left hemisphere. (b) Right hemisphere.

Figure 28:

Electrode locations in patient 3. (a) Left hemisphere. (b) Right hemisphere.

Plots of MI, CC, EC, and latent inputs measures from the SLRPM method for seven seizures of this patient are shown in Figures 29, 30, 31, and 32, respectively. Ictal dynamics is clearly visible from these plots and similar to that in patients 1 and 2.

Figure 29:

MI plots for the SLRPM method for seizures 1, 2, 3, 4, 5, 6, and 7, shown in panels a, b, c, d, e, f, and g, respectively, in patient 3. Green lines in all plots denote the seizure onset time.

Figure 29:

MI plots for the SLRPM method for seizures 1, 2, 3, 4, 5, 6, and 7, shown in panels a, b, c, d, e, f, and g, respectively, in patient 3. Green lines in all plots denote the seizure onset time.

Figure 30:

CC plots for the SLRPM method for seizures 1, 2, 3, 4, 5, 6, and 7, shown in panels a, b, c, d, e, f, and g, respectively, in patient 3. Green lines in all plots denote the seizure onset time.

Figure 30:

CC plots for the SLRPM method for seizures 1, 2, 3, 4, 5, 6, and 7, shown in panels a, b, c, d, e, f, and g, respectively, in patient 3. Green lines in all plots denote the seizure onset time.

Figure 31:

EC plots for the SLRPM method for seizures 1, 2, 3, 4, 5, 6, and 7, shown in panels a, b, c, d, e, f, and g, respectively, in patient 3. Green lines in all plots denote the seizure onset time.

Figure 31:

EC plots for the SLRPM method for seizures 1, 2, 3, 4, 5, 6, and 7, shown in panels a, b, c, d, e, f, and g, respectively, in patient 3. Green lines in all plots denote the seizure onset time.

Figure 32:

Plots of the latent inputs estimated by the SLRPM method for seizures 1, 2, 3, 4, 5, 6, and 7, shown in panels a, b, c, d, e, f, and g, respectively, in patient 3. Green lines in all plots denote the seizure onset time.

Figure 32:

Plots of the latent inputs estimated by the SLRPM method for seizures 1, 2, 3, 4, 5, 6, and 7, shown in panels a, b, c, d, e, f, and g, respectively, in patient 3. Green lines in all plots denote the seizure onset time.

### 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.

Figure 33:

Electrode locations in patient 4.

Figure 33:

Electrode locations in patient 4.

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.

Figure 34:

MI plots for the SLRPM method for seizures 1, 2, and 3, shown in panels a, b, and c, respectively, in patient 4. Green lines in all plots denote the seizure onset time.

Figure 34:

MI plots for the SLRPM method for seizures 1, 2, and 3, shown in panels a, b, and c, respectively, in patient 4. Green lines in all plots denote the seizure onset time.

Figure 35:

CC plots for the SLRPM method for seizures 1, 2, and 3, shown in panels a, b, and c, respectively, in patient 4. Green lines in all plots denote the seizure onset time.

Figure 35:

CC plots for the SLRPM method for seizures 1, 2, and 3, shown in panels a, b, and c, respectively, in patient 4. Green lines in all plots denote the seizure onset time.

Figure 36:

EC plots for the SLRPM method for seizures 1, 2, and 3, shown in panels a, b, and c, respectively, in patient 4. Green lines in all plots denote the seizure onset time.

Figure 36:

EC plots for the SLRPM method for seizures 1, 2, and 3, shown in panels a, b, and c, respectively, in patient 4. Green lines in all plots denote the seizure onset time.

Figure 37:

Plots of the latent inputs estimated by the SLRPM method for seizures 1, 2, and 3, shown in panels a, b, and c, respectively, in patient 4. Green lines in all plots denote the seizure onset time.

Figure 37:

Plots of the latent inputs estimated by the SLRPM method for seizures 1, 2, and 3, shown in panels a, b, and c, respectively, in patient 4. Green lines in all plots denote the seizure onset time.

### 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.

Figure 38:

Electrode locations in patient 5.

Figure 38:

Electrode locations in patient 5.

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.

Figure 39:

MI plots for the SLRPM method for four seizures in patient 5. Green lines in all plots denote the seizure onset time.

Figure 39:

MI plots for the SLRPM method for four seizures in patient 5. Green lines in all plots denote the seizure onset time.

Figure 40:

CC plots for the SLRPM method for four seizures in patient 5. Green lines in all plots denote the seizure onset time.

Figure 40:

CC plots for the SLRPM method for four seizures in patient 5. Green lines in all plots denote the seizure onset time.

Figure 41:

EC plots for the SLRPM method for four seizures in patient 5. Green lines in all plots denote the seizure onset time.

Figure 41:

EC plots for the SLRPM method for four seizures in patient 5. Green lines in all plots denote the seizure onset time.

Figure 42:

Plots of the latent inputs estimated by the SLRPM method for four seizures in patient 5. Green lines in all plots denote the seizure onset time.

Figure 42:

Plots of the latent inputs estimated by the SLRPM method for four seizures in patient 5. Green lines in all plots denote the seizure onset time.

## 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.

## References

Aguirre
,
A. G.
,
Portes
,
L. L.
, &
Letellier
,
C.
(
2017
).
Observability and synchronization of neuron models
.
Chaos
,
27
(
10
),
103103
.
Allen
,
E. A.
,
Damaraju
,
E.
,
Plis
,
S. M.
,
Erhardt
,
E. B.
,
Eichele
,
E.
, &
Calhoun
,
V. D.
(
2012
).
Tracking whole-brain connectivity dynamics in the resting state
.
Cerebral Cortex
,
24
,
663
676
.
Anand
,
A.
,
Li
,
Y.
,
Wang
,
Y.
,
Wu
,
J.
,
Gao
,
S.
,
Bukhari
,
L.
, …
Lowe
,
M. J.
(
2005
).
Antidepressant effect on connectivity of the mood-regulating circuit: An fMRI study
.
Neuropsychopharmacology
,
30
(
7
),
1334
1344
.
Attwell
,
D.
, &
,
C.
(
2002
).
The neural basis of functional brain imaging signals
.
Trends in Neurosciences
,
25
(
12
),
621
625
.
Banerjee
,
O.
,
El Ghaoui
,
L.
, &
d'Aspremont
,
A.
(
2008
).
Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data
.
Journal of Machine Learning Research
,
9
,
485
516
.
Bassett
,
D. S.
,
Greenfield
,
D. L.
,
Lindenberg
,
A. M.
,
Weinberger
,
D. R.
,
Moore
,
S. W.
, &
Bullmore
,
E. T.
(
2010
).
Efficient physical embedding of topologically complex information processing networks in brains and computer circuits
.
PLoS Computational Biology
,
6
(
4
),
e1000748
.
Bassett
,
D. S.
,
Wymbs
,
N. F.
,
Porter
,
M. A.
,
Mucha
,
P. J.
,
Carlson
,
J. M.
, &
Grafton
,
S. T.
(
2011
).
Dynamic reconfiguration of human brain networks during learning
.
Proceedings of the National Academy of Sciences USA
,
108
(
18
),
7641
7646
.
Biswal
,
B.
,
Yetkin
,
F. Z.
,
Haughton
,
V. M.
, &
Hyde
,
J. S.
(
1995
).
Functional connectivity in the motor cortex of resting human brain using echo-planar MRI
.
Magnetic Resonance in Medicine
,
34
(
4
),
537
541
.
Blondel
,
V. D.
,
Guillaume
,
J. L.
,
Lambiotte
,
R.
, &
Lefebvre
,
E.
(
2008
).
Fast unfolding of communities in large networks
.
Journal of Statistical Mechanics: Theory and Experiment
,
2008
(
10
),
P10008
.
Blumenfeld
,
H.
(
2012
).
Impaired consciousness in epilepsy
.
Lancet Neurology
,
11
(
9
),
814
826
.
Brandes
,
U.
,
Delling
,
D.
,
Gaertler
,
M.
,
Gorke
,
R.
,
Hoefer
,
M.
,
Nikoloski
,
Z.
, &
Wagner
,
D.
(
2008
).
On modularity clustering
.
IEEE Transactions on Knowledge and Data Engineering
,
20
(
2
),
172
188
.
Bruno
,
J. L.
,
Hosseini
,
S. M. H.
,
Saggar
,
M.
,
Quintin
,
E.
,
Raman
,
M. M.
, &
Reiss
,
A. L.
(
2017
).
Altered brain network segregation in fragile X syndrome revealed by structural connectomics
.
Cerebral Cortex
,
27
(
3
),
2249
2259
.
Burns
,
S. P.
,
Santaniello
,
S.
,
Yaffe
,
R. B.
,
Jouny
,
C. C.
,
Crone
,
N. E.
,
Bergey
,
G. K.
, …
Sarma
,
S. V.
(
2014
).
Network dynamics of the brain and influence of the epileptic seizure onset zone
.
Proceedings of the National Academy of Sciences USA
,
111
(
49
),
E5321
E5330
.
Chandrasekaran
,
V.
,
Parrilo
,
P. A.
, &
Willsky
,
A. S.
(
2012
).
Latent variable graphical model selection via convex optimization
.
Annals of Statistics
,
40
(
4
),
1935
1967
.
Chu
,
C. J.
,
Tanaka
,
N.
,
Diaz
,
J.
,
Edlow
,
B. L.
,
Wu
,
O.
,
Hamalainen
, M.,
… Kramer
, M. A.
(
2015
).
EEG functional connectivity is partially predicted by underlying white matter connectivity
.
NeuroImage
,
108
,
23
33
.
Cimenser
,
A.
,
Purdon
,
P. L.
,
Pierce
,
E. T.
,
Walsh
,
J. L.
,
Salazar-Gomez
,
A. F.
,
Harrell
,
P. G.
, …
Brown
,
E. N.
(
2011
).
Tracking brain states under general anesthesia by using global coherence analysis
.
Proceedings of the National Academy of Sciences USA
,
108
(
21
),
8832
8837
.
Cole
,
M. W.
,
Bassett
,
D. S.
,
Power
,
J. D.
,
Braver
,
T. S.
, &
Petersen
,
S. E.
(
2014
).
Intrinsic and task-evoked network architectures of the human brain
.
Neuron
,
83
(
1
),
238
251
.
Cordes
,
D.
,
Haughton
,
V. M.
,
Arfanakis
,
K.
,
Wendt
,
G. J.
,
Turski
,
P. A.
,
Moritz
,
C. H.
, …
Meyerand
,
M. E.
(
2000
).
Mapping functionally related regions of brain with functional connectivity MR imaging
.
,
21
(
9
),
1636
1644
.
Das
,
A.
(
2018
).
The importance of latent inputs for analyzing the human brain
.
PhD diss., University of California, San Diego
.
Das
,
A.
,
Cash
,
S. S.
, &
Sejnowski
,
T. J.
(
2018
).
Heterogeneity of human epileptic seizures
.
Manuscript in preparation
.
Das
,
A.
,
Sampson
,
A. L.
,
Lainscsek
,
C.
,
Muller
,
L.
,
Lin
,
W.
,
Doyle
,
J. C.
, …
Sejnowski
,
T. J.
(
2017
).
Interpretation of the precision matrix and its application in estimating sparse brain connectivity during sleep spindles from human electrocorticography recordings
.
Neural Computation
,
29
(
3
),
603
642
.
Dempster
,
A. P.
(
1972
).
Covariance selection
.
Biometrics
,
28
(
1
),
157
175
.
Diessen
,
E. V.
,
Diederen
,
S. J. H.
,
Braun
,
K. P. J.
,
Jansen
,
F. E.
, &
Stam
,
C. J.
(
2013
).
Functional and structural brain networks in epilepsy: What have we learned?
Epilepsia
,
54
(
11
),
1855
1865
.
Fox
,
M. D.
,
Corbetta
,
M.
,
Snyder
,
A. Z.
,
Vincent
,
J. L.
, &
Raichle
,
M. E.
(
2006
).
Spontaneous neuronal activity distinguishes human dorsal and ventral attention systems
.
Proceedings of the National Academy of Sciences USA
,
26
(
103
),
10046
10051
.
Fox
,
M. D.
,
Snyder
,
A. Z.
,
Vincent
,
J. L.
,
Corbetta
,
M.
,
Van Essen
,
D. C.
, &
Raichle
,
M. E.
(
2005
).
The human brain is intrinsically organized into dynamic, anticorrelated functional networks
.
Proceedings of the National Academy of Sciences USA
,
102
(
27
),
9673
9678
.
Friedman
,
J.
,
Hastie
,
T.
, &
Tibshirani
,
R.
(
2008
).
Sparse inverse covariance estimation with the graphical lasso
.
Biostatistics
,
9
(
3
),
432
441
.
Giorge
,
F. S.
,
Pizzanelli
,
C.
,
Biagioni
,
F.
,
Murri
,
L.
, &
Fornia
,
F.
(
2004
).
The role of norepinephrine in epilepsy: From the bench to the bedside
.
Neuroscience and Biobehavioral Reviews
,
28
(
5
),
507
524
.
Glasser
,
M. F.
,
Coalson
,
T. S.
,
Robinson
,
E. C.
,
Hacker
,
C. D.
Harwell
,
J.
,
Yacoub
,
E.
, …
Van Essen
,
D. C.
(
2016
).
A multi-modal parcellation of human cerebral cortex
.
Nature
,
536
,
171
178
.
Greicius
,
M. D.
,
Krasnow
,
B.
,
Reiss
,
A. L.
, &
Menon
,
V.
(
2003
).
Functional connectivity in the resting brain: A network analysis of the default mode hypothesis
.
Proceedings of the National Academy of Sciences USA
,
1
(
100
),
253
258
.
Hampson
,
M.
,
Peterson
,
B. S.
,
Skudlarski
,
P.
,
Gatenby
,
J. C.
, &
Gore
,
J. C.
(
2002
).
Detection of functional connectivity using temporal correlations in MR images
.
Human Brain Mapping
,
15
(
4
),
247
262
.
Honey
,
C. J.
,
Sporns
,
O.
,
Cammoun
,
L.
,
Gigandet
,
X.
,
Thiran
,
J. P.
,
Meuli
,
R.
, &
Hagmann
,
P.
(
2009
).
Predicting human resting-state functional connectivity from structural connectivity
.
Proceedings of the National Academy of Sciences USA
,
106
(
6
),
2035
2040
.
Hsieh
,
C. J.
,
Sustik
,
M. A.
,
Dhillon
,
I. S.
, &
Ravikumar
,
P.
(
2011
). Sparse inverse covariance matrix estimation using quadratic approximation. In
L.
Shawe-Taylor
,
R. S.
Zemel
,
P. L.
Bartlett
,
F.
Pereira
, &
K. Q.
Weinberger
(Eds.),
Advances in neural information processing systems
,
24
(pp.
2330
2338
).
Red Hook, NY
:
Curran
.
Hsieh
,
C. J.
,
Sustik
,
M. A.
,
Dhillon
,
I. S.
,
Ravikumar
,
P.
, &
Poldrack
,
R.
(
2013
). BIG & QUIC: Sparse inverse covariance estimation for a million variables. In
C. J. C.
Burges
,
L.
Bottou
,
M.
Welling
,
Z.
Ghahramani
, &
K. Q.
Weinberger
(Eds.),
Advances in neural information processing systems
,
26
(pp.
3165
3173
).
Red Hook, NY
:
Curran
.
Huth
,
A. G.
,
de Heer
,
W. A.
,
Griffiths
,
T. A.
,
Theunissen
,
F. E.
, &
Gallant
,
J. L.
(
2016
).
Natural speech reveals the semantic maps that tile human cerebral cortex
.
Nature
,
532
,
453
458
.
Jerger
,
K. K.
,
Netoff
,
T. I.
,
Francis
,
J. T.
,
Sauer
,
T.
,
Pecora
,
L.
,
Weinstein
,
S. L.
, &
Schiff
,
S. J.
(
2001
).
Early seizure detection
.
Journal of Clinical Neurophysiology
,
18
(
3
),
259
268
.
Joyce
,
K. E.
,
Laurienti
,
P. J.
,
Burdette
,
J. H.
, &
Hayasaka
,
S.
(
2010
).
A new measure of centrality for brain networks
.
PloS One
,
5
(
8
),
e12200
.
Kramer
,
M. A.
, &
Cash
,
S. S.
(
2012
).
Epilepsy as a disorder of cortical network organization
.
Neuroscientist
,
18
(
4
),
360
372
.
Kramer
,
M. A.
,
Eden
,
U. T.
,
Kolaczyk
,
E. D.
,
Zepeda
,
R.
,
Eskandar
,
E. N.
, &
Cash
,
S. S.
(
2010
).
Coalescence and fragmentation of cortical networks during focal seizures
.
Journal of Neuroscience
,
30
(
30
),
10076
10085
.
Krimer
,
L. S.
,
Muly
,
E. C. III
,
Williams
,
G. V.
, &
Goldman-Rakic
,
P. S.
(
1998
).
Dopaminergic regulation of cerebral cortical microcirculation
.
Nature Neuroscience
,
1
(
4
),
286
289
.
,
F. A.
, &
Moshe, S.
L.
(
2008
).
How do seizures stop?
Epilepsia
,
49
(
10
),
1651
1664
.
Laughlin
,
S. B.
, &
Sejnowski
,
T. J.
(
2003
).
Communication in neuronal networks
.
Science
,
301
(
5641
),
1870
1874
.
Lauritzen
,
S. L.
(
1996
).
Graphical models
.
New York
:
Oxford University Press
.
Liu
,
Y. Y.
,
Slotine
,
J. J.
, &
Barabasi
,
A. L.
(
2013
).
Observability of complex systems
.
Proceedings of the National Academy of Sciences USA
,
110
(
7
),
2460
2465
.
Lohmann
,
G.
,
Margulies
,
D. S.
,
Horstmann
,
A.
,
Pleger
,
B.
,
Lepsien
,
J.
,
Goldhahn
,
D.
, …
Turner
,
R.
(
2010
).
Eigenvector centrality mapping for analyzing connectivity patterns in fMRI data of the human brain
.
PloS One
,
5
(
4
),
e10232
.
Lynall
,
M.
,
Bassett
,
D. S.
,
Kerwin
,
R.
,
McKenna
,
P. J.
,
Kitzbichler
,
M.
,
Muller
,
U.
, &
Bullmore
,
E.
(
2010
).
Functional connectivity and brain networks in schizophrenia
.
Journal of Neuroscience
,
30
(
28
),
9477
9487
.
Ma
,
S.
,
Xue
,
L.
, &
Zou
,
H.
(
2013
).
Alternating direction methods for latent variable gaussian graphical model selection
.
Neural Computation
,
25
(
8
),
2172
2198
.
Marrelec
,
G.
,
Horwitz
,
B.
,
Kim
,
J.
,
Pelegrini-Issac
,
M.
,
Benali
,
H.
, &
Doyon
,
J.
(
2007
).
Using partial correlation to enhance structural equation modeling of functional MRI data
.
Magnetic Resonance Imaging
,
25
(
8
),
1181
1189
.
Meunier
,
D.
,
Lambiotte
,
R.
, &
Bullmore
,
E. T.
(
2010
).
Modular and hierarchically modular organization of brain networks
.
Frontiers in Neuroscience
,
4
(
200
).
Meyer-Base
,
A.
,
Roberts
,
R. G.
,
Illan
,
I. A.
,
Meyer-Base
,
U.
,
Lobbes
,
M.
,
,
A.
, &
Pinker-Domenig
,
A.
(
2017
).
Dynamical graph theory networks methods for the analysis of sparse functional connectivity networks and for determining pinning observability in brain networks
.
Frontiers in Computational Neuroscience
,
11
(
87
).
Monti
,
R. P.
,
Hellyer
,
P.
,
Sharp
,
D.
,
Leech
,
R.
,
Anagnostopoulos
,
C.
, &
Montana
,
G.
(
2014
).
Estimating time-varying brain connectivity networks from functional MRI time series
.
NeuroImage
,
103
,
427
443
.
Newman
,
M. E. J.
(
2006
).
Modularity and community structure in networks
.
Proceedings of the National Academy of Sciences USA
,
103
(
23
),
8577
8582
.
Newman
,
M. E. J.
(
2010
).
Networks: An introduction
.
Oxford
:
Oxford University Press
.
Newman
,
M. E. J.
, &
Girvan
,
M.
(
2004
).
Finding and evaluating community structure in networks
.
Physical Review E
,
69
(
2
),
026113
.
Nunez
,
P. L.
, &
Srinivasan
,
R.
(
2006
).
Electric fields of the brain: The neurophysics of EEG
.
New York
:
Oxford University Press
.
Onnela
,
J. P.
,
Saramaki
,
J.
,
Kertesz
,
J.
, &
,
K.
(
2005
).
Intensity and coherence of motifs in weighted complex networks
.
Physical Review E
,
71
(
6
),
764
772
.
Oztoprak
,
F.
,
Nocedal
,
J.
,
Rennie
,
S.
, &
Olsen
,
P. A.
(
2012
). Newton-like methods for sparse inverse covariance estimation. In
F.
Pereira
,
C. J. C.
Burges
,
L.
Bottou
, &
K. Q.
Weinberger
(Eds.),
Advances in neural information processing systems, 25
(pp.
764
772
).
Red Hook, NY
:
Curran
.
Ponten
,
S. C.
,
Bartolomei
,
F.
, &
Stam
,
C. J.
(
2007
).
Small-world networks and epilepsy: Graph theoretical analysis of intracerebrally recorded mesial temporal lobe seizures
.
Clinical Neurophysiology
,
118
(
4
),
918
927
.
Reichardt
,
J.
, &
Bornholdt
,
S.
(
2006
).
Statistical mechanics of community detection
.
Physical Review E
,
74
(
1
),
016110
.
Ronhovde
,
P.
, &
Nussinov
,
Z.
(
2009
).
Multiresolution community detection for megascale networks by information-based replica correlations
.
Physical Review E
,
80
(
1
),
016109
.
Rosa
,
M. J.
,
Portugal
,
L.
,
Hahn
,
T.
,
Fallgatter
,
A. J.
,
Garrido
,
M. I.
,
Taylor
,
J. S.
, &
Mourao-Miranda
,
J.
(
2015
).
Sparse network-based models for patient classification using fMRI
.
NeuroImage
,
105
,
493
506
.
Rothman
,
A. J.
,
Bickel
,
P. J.
,
Levina
,
E.
, &
Zhu
,
J.
(
2008
).
Sparse permutation invariant covariance estimation
.
Electronic Journal of Statistics
,
2
,
494
515
.
Rubinov
,
M.
, &
Sporns
,
O.
(
2010
).
Complex network measures of brain connectivity: Uses and interpretations
.
NeuroImage
,
52
(
3
),
1059
1069
.
Rubinov
,
M.
, &
Sporns
,
O.
(
2011
).
Weight-conserving characterization of complex functional brain networks
.
NeuroImage
,
56
(
4
),
2068
2079
.
Ryali
,
S.
,
Chen
,
T.
,
Supekar
,
K.
, &
Menon
,
V.
(
2012
).
Estimation of functional connectivity in fMRI data using stability selection-based sparse partial correlation with elastic net penalty
.
NeuroImage
,
59
(
4
),
3852
3861
.
Saggar
,
M.
,
Hosseini
,
S. M. H.
,
Bruno
,
J. L.
,
Quintin
,
E.
,
Raman
,
M. M.
,
Kesler
,
S. R.
, &
Reiss
,
A. L.
(
2015
).
Estimating individual contribution from group-based structural correlation networks
.
NeuroImage
,
120
,
274
184
.
,
R.
,
Suckling
,
J.
,
Coleman
,
M. R.
,
Pickard
,
J. D.
,
Menon
,
D.
, &
Bullmore
,
E.
(
2005
).
Neurophysiological architecture of functional magnetic resonance images of human brain
.
Cerebral Cortex
,
15
,
1332
1342
.
Saramaki
,
J.
,
Kivela
,
M.
,
Onnela
,
J. P.
,
,
K.
, &
Kertesz
,
J.
(
2007
).
Generalizations of the clustering coefficient to weighted complex networks
.
Physical Review E
,
75
(
2
),
027105
.
Sato
,
A.
, &
Sato
,
Y.
(
1992
).
Regulation of regional cerebral blood flow by cholinergic fibers originating in the basal forebrain
.
Neuroscience Research
,
14
(
4
),
242
274
.
Schindler
,
K. A.
,
Bialonski
,
S.
,
Horstmann
,
M. T.
,
Elger
,
C. E.
, &
Lehnertz
,
K.
(
2008
).
Evolving functional network properties and synchronizability during human epileptic seizures
.
Chaos
,
18
(
3
),
033119
.
Scheinberg
,
K.
,
Ma
,
S.
, &
Goldfarb
,
D.
(
2010
). Sparse inverse covariance selection via alternating linearization methods. In
J. D.
Lafferty
,
C. K. I.
Williams
,
J.
Shawe-Taylor
,
R. S.
Zemel
, &
A.
Culotta
(Eds.),
Advances in neural information processing systems
,
23
(pp.
2101
2109
).
Red Hook, NY
:
Curran
.
Shafi
,
M. M.
,
Westover
,
M. B.
,
Oberman
,
L.
,
Cash
,
S. S.
, &
Pascual-Leone
,
A.
(
2014
).
Modulation of EEG functional connectivity networks in subjects undergoing repetitive transcranial magnetic stimulation
.
Brain Topography
,
27
(
1
),
172
191
.
Siegle
,
G. J.
,
Thompson
,
W.
,
Carter
,
C. S.
,
Steinhauer
,
S. R.
, &
Thase
,
M. E.
(
2007
).
Increased amygdala and decreased dorsolateral prefrontal BOLD responses in unipolar depression: Related and independent features
.
Biological Psychiatry
,
61
(
2
),
198
209
.
Smith
,
S. M.
,
Miller
,
K. L.
,
Salimi-Khorshidi
,
G.
,
Webster
,
M.
,
Beckmann
,
C. F.
,
Nichols
,
T. E.
, …
Woolrich
,
M. W.
(
2011
).
Network modelling methods for fMRI
.
NeuroImage
,
54
(
2
),
875
891
.
Sporns
,
O.
(
2010
).
Brain connectivity toolbox
Sporns
,
O.
(
2011
).
The nonrandom brain: Efficiency, economy, and complex dynamics
.
Frontiers in Computational Neuroscience
,
5
(
5
).
Stam
,
C. J.
(
2005
).
Nonlinear dynamical analysis of EEG and MEG: Review of an emerging field
.
Clinical Neurophysiology
,
116
(
10
),
2266
2301
.
Stam
,
C. J.
(
2014
).
Modern network science of neurological disorders
.
Nature Reviews Neuroscience
,
15
,
683
695
.
,
M.
,
McCormick
,
D. A.
, &
Sejnowski
,
T. J.
(
1993
).
Thalamocortical oscillations in the sleeping and aroused brain
.
Science
,
262
(
5134
),
679
685
.
Sun
,
Y.
,
Danila
,
B.
,
Josić
,
K.
, &
Bassler
,
K. E.
(
2009
).
Improved community structure detection using a modified fine-tuning strategy
.
Europhysics Letters
,
86
(
2
),
28004
.
Supekar
,
K.
,
Menon
,
V.
,
Rubin
,
D.
,
Muse
,
M.
, &
Greicius
,
M. D.
(
2008
).
Network analysis of intrinsic functional brain connectivity in Alzheimer's disease
.
PLoS Computational Biology
,
4
(
6
),
e1000100
.
Uddin
,
L. Q.
,
Kelly
,
A. C.
,
Biswal
,
B. B.
,
Castellanos
,
F. X.
, &
Milham
,
M. P.
(
2009
).
Functional connectivity of default mode network components: Correlation, anticorrelation, and causality
.
Human Brain Mapping
,
30
(
2
),
625
637
.
Varoquaux
,
G.
,
Gramfort
,
A.
,
Jean-Baptiste
,
P.
, &
Thirion
,
B.
(
2010
). Brain covariance selection: Better individual functional connectivity models using population prior. In
J. D.
Lafferty
,
C. K. I.
Williams
,
J.
Shawe-Taylor
,
R. S.
Zemel
, &
A.
Culotta
(Eds.),
Advances in neural information processing systems
,
23
(pp.
2334
2342
).
Red Hook, NY
:
Curran
.
Varsavsky
,
A.
,
Mareels
,
I.
, &
Cook
,
M.
(
2010
).
Epileptic seizures and the EEG: Measurement, models, detection and prediction
.
Boca Raton, FL
:
CRC Press
.
Vega-Zelaya
,
L.
,
Pastor
,
J.
,
de Sola
,
R. G.
, &
Ortega
,
G. J.
(
2015
).
Disrupted ipsilateral network connectivity in temporal lobe epilepsy
.
PloS One
,
10
(
10
),
e0140859
.
Vertes
,
P. E.
,
Alexander-Bloch
,
A. F.
,
Gogtay
,
N.
,
Giedd
,
J. N.
,
Rapoport
,
J. L.
, &
Bullmore
,
E. T.
(
2012
).
Simple models of human brain functional networks
.
Proceedings of the National Academy of Sciences
,
109
(
15
),
5868
5873
.
Wang
,
Y.
,
Ghumare
,
E.
,
Vandenberghe
,
R.
, &
Dupont
,
P.
(
2017
).
Comparison of different generalizations of clustering coefficient and local efficiency for weighted undirected graphs
.
Neural Computation
,
29
(
2
),
313
331
.
Wang
,
Y.
,
Kang
,
J.
,
Kemmer
,
P. B.
, &
Guo
,
Y.
(
2016
).
An efficient and reliable statistical method for estimating functional connectivity in large scale brain networks using partial correlation
.
Frontiers in Neuroscience
,
10
(
123
).
Whalen
,
A. J.
,
Brennan
,
S. N.
,
Sauer
,
T. D.
, &
Schiff
,
S. J.
(
2015
).
Observability and controllability of nonlinear networks: The role of symmetry
.
Physical Review X
,
5
(
1
),
011005
.
Whittaker
,
J.
(
1990
).
Graphical models in applied multivariate statistics
.
Hoboken, NJ
:
Wiley
.
Xiong
,
J.
,
Parsons
,
L. M.
,
Gao
,
J. H.
, &
Fox
,
P. T.
(
1999
).
Interregional connectivity to primary motor cortex revealed using MRI resting state images
.
Human Brain Mapping
,
8
(
2–3
),
151
156
.
Yatsenko
,
D.
,
Josić
,
K.
,
Ecker
,
A. S.
,
Froudarakis
,
E.
,
Cotton
,
R. J.
, &
Tolias
,
A. S.
(
2015
).
Improved estimation and interpretation of correlations in neural circuits
.
PLoS Computational Biology
,
11
(
3
),
e1004083
.
Yuan
,
M.
, &
Lin
,
Y.
(
2007
).
Model selection and estimation in the gaussian graphical model
.
Biometrika
,
94
(
1
),
19
35
.
Zhou
,
D.
,
Thompson
,
W. K.
, &
Siegle
,
G.
(
2009
).
MATLAB toolbox for functional connectivity
.
NeuroImage
,
47
(
4
),
1590
1607
.
Zuo
,
X. N.
,
Ehmke
,
R.
,
Mennes
,
M.
,
Imperati
,
D.
,
Castellanos
,
F. X.
,
Sporns
,
O.
, &
Milham
,
M. P.
(
2012
).
Network centrality in the human functional connectome
.
Cerebral Cortex
,
22
(
8
),
1862
1875
.