Modeling Brain Connectivity with Graphical Models on Frequency Domain

Modeling Brain Connectivity with Graphical Models on Frequency Domain

Xu Gao Department of Statistics, University of California, Irvine, CA, USA Weining Shen Department of Statistics, University of California, Irvine, CA, USA Chee-Ming Ting Steven C. Cramer Department of Anatomy and Neurobiology, University of California, Irvine, CA, USA Ramesh Srinivasan Department of Cognitive Sciences, University of California, Irvine, CA, USA Hernando Ombao Statistics Program, King Abdullah University of Science and Technology, Thuwal, Saudi Arabia

Multichannel electroencephalograms (EEGs) have been widely used to study cortical connectivity during acquisition of motor skills. In this paper, we introduce copula Gaussian graphical models on spectral domain to characterize dependence in oscillatory activity between channels. To obtain a simple and robust representation of brain connectivity that can explain the most variation in the observed signals, we propose a framework based on maximizing penalized likelihood with Lasso regularization to search for the sparse precision matrix. To address the optimization problem, graphical Lasso, Ledoit-Wolf and sparse estimation of a covariance matrix (SPCOV) algorithms were modified and implemented. Simulations show the benefit of using the proposed algorithms in terms of robustness and small estimation errors. Furthermore, analysis of the EEG data in a motor skill task conducted using algorithms of modified graphical LASSO and Ledoit-Wolf, reveal a sparse pattern of brain connectivity among cortices which is consistent with the results from other work in the literature.

Keywords:  Brain connectivity; Copula Gaussian graphical model; graphical LASSO; SPCOV; Ledoit-Wolf algorithm; electroencephalograph.

1 Introduction

A graph is a model representation of a complex system determined by a set of nodes (vertices) and edges connecting them (Trudeau, 1993). On the foundation of graph and probability theory, graphical models (probabilistic graphical model) are widely used in the communities of Bayesian statistics, statistical learning and machine learning. In the framework of graphical models, each node represents a random variable and edges denote the probabilistic relationship between nodes. The graph depicts the structure where the joint distribution of random variables can be decomposed into a product of factors depending only on subsets of variables (Bishop, 2006).

Graphs have been introduced to model brain connectivity through the way that nodes represent cortical and subcortical regions while edges can be denoted as functional and structural connections between cortical nodes (Bullmore& Bassett, 2011). In the literature of brain graph models, much work has been done for all major modalities of functional magnetic resonance imaging (fMRI) and neurophysiological data. To name a few, functional brain graphs and their relevant works have been constructed from fMRI (Achard& Bullmore, 2007), electroencephalography (EEG) signals (Michelovannis et al., 2006; Kemmer et al., 2015), magnetoencephalography (MEG) data (Stam, 2004) and Local Field Potentials (LFPs) signals (Gao et al., 2016, 2018). From diffusion tensor imaging (DTI) and diffusion spectrum imaging(DSI), structural brain graphs have been studied by Gong et al. (2008). Among them, sparse graphical models, which are widely discussed in Friedman et al. (2008), are highly efficient in inferring multielectrode brain recordings. The sparsity of graph provides a robust approach that highlights the most significant interactions between brain cortices and helps to interpret the data (Dauwels et al., 2012A).

Motivated by the advantages of using sparse graph, we considered the problem of modeling brain connectivity from EEG data under the framework of sparse graphical models. Inspired by the work of Dauwels et al. (2012B), the main contributions of the paper are as follows: (1.) We introduce copula Gaussian graphical models to account for the non-Gaussianity of signals on frequency domain; (2.) We develop a framework to capture the between-channel dependence in the oscillatory activity; (3.) By including a regularization term, we are able to obtain a simple and robust representation that uncovers the most critical connectivities between brain cortices.

In this paper, the core methodology boils down to the sparse estimation of covariance matrices. Specifically, we regularize the log likelihood with a lasso penalty on the entries of covariance matrix as the objective function. The penalty is used to reduce the effective number of brain connectivity and thus produces sparse and robust estimates (Bein& Tibshirani, 2010). To solve the optimization problem, many algorithms have been introduced in the literature. Rothman et al. (2008) propose a novel algorithm with the assumption of ordering to the variables. Butte et al. (2000) introduce relevance work in working on the optimization problem in the way that pairwise correlation beyond a threshold are linked by an edge. Rothman et al. (2009) propose an algorithm by introducing shrinkage operators. Rothman et al. (2010) utilize lasso-regression based method to solve the optimization problem. Throughout this paper, we develop modified algorithms based on the works of graphical lasso (Friedman et al., 2008), sparse estimation of covariance proposed by Bein& Tibshirani (2010) and LedoitWolf algorithm (Ledoit& Wolf, 2012).

The paper is organized as follows. In section 2, we briefly review copula Gaussian graphical models and its application in EEG data. In section 3, we formulate the optimization problem and propose three algorithms to solve it. In section 4, we conduct simulations to test the robustness and performance of the three algorithms. In section 4, we apply the EEG data to the optimization problem in searching for the sparse connectivity matrix.

2 Background on Graphical Models in Brain Connectivity

In this section, we discuss the preliminaries on graphical models and its application in modeling brain connectivity from EEG data.

2.1 Copula Gaussian Graphical Model

Suppose we have non-Gaussian random variables . We define hidden Gaussian random variables through the relationship that (Dauwels et al., 2012A)


where is the precision matrix, is the cumulative distribution function of Gaussian random variables and is the empirical cumulative distribution function of . In practice, can be estimated by

2.2 Graphical Models in EEG data

In fact, the EEG data obtained from different electrodes is multiple non-Gaussian time series. By implementing Copula Gaussian Graphical Models, the original time series have been transformed into Gaussian dataset. In the following work, we implement methodology in graphical models for multivariate Gaussian time series.

Suppose a graphical model uniquely determines the conditional independence on Gaussian process In graph , each node denotes a single time series The absence of edge between and denotes the conditional independence between time series and . Under the assumption that the cross-variance function of is summable,

we define the cross spectral density matrix of as

where denotes the Fourier transform. As a result of Dahlhaus (2000), the Gaussian process and are conditional independence if and only if

In practice, we use as the empirical variance-covariance matrix of the time series

3 Optimization Problem Statement and the Proposed Algorithms

Following the methodology discussed in section 2, we have transformed the EEG data into quasi-Gaussian time series with empirical variance-covariance matrix In this section, we formulate the final optimization problem in modeling the brain connectivity and propose three algorithms in solving the problem.

3.1 Optimization Problem Formulation

In sparse graphical models, true brain connectivity involving the strongest and the most relevant connections is uniquely determined by the sparse precision matrix (the inverse of covariance matrix). The objective function, defined as the regularized negative log-likelihood function is given by

subject to

where is the empirical variance-covariance matrix discussed in section 2, is the norm which is given by the sum over the absolute values of entries in matrix and is a tuning parameter controlling the amount of shrinkage.

3.2 The Proposed Algorithms

We propose three algorithms to address the optimization problem (A-3).

SPCOV (Majorize-Minimize) algorithm

In the objective function (A-3), is convex while is concave, thus a majorize-minimize scheme could be used. In summary, this algorithm presents as two loops, the outer loop approximates the non-convex problem and the inner loop solves each convex relaxation.

Graphical lasso algorithm

We also propose a modified graphical lasso algorithm. The rationale is that suppose is the estimate of , then one can solve the problem by optimizing over each row and corresponding column of in a block coordinate descent fashion.

Modified Ledoit-Wolf algorithm

Inspired by the works of Ledoit& Wolf (2012), Fiecas et al. (2010) and Fiecas& Ombao (2011), we also implement a modified Ledoit-Wolf algorithm to address the optimization problem A-3. Specifically, we utilize the sample covariance and the maximum likelihood estimator to obtain a shrinkage estimator which compromises between variance and bias. The shrinkage intensity is evaluated by minimizing a risk function based on mean square errors.

4 Simulations

In this section, two simulation scenarios were considered to evaluate the performance of the proposed algorithms. In each scenario, we created different types of sparse symmetric positive definite matrices to be the true covariance structure. We randomly generated 200 samples from the true covariance structure and then used the data as the input of the proposed algorithms. Each scenarios was repeated 1000 times. We evaluated the results by its ability to correctly identify the zero elements of and the discrepancy with the true covariance matrix by root- mean-square error, and entropy loss, respectively, where is the number of parameters.

The first scenario was from Cliques Model. In this model, we set to be the covariance matrix. Each in the diagonal represented a dense matrix. Other parts of the matrix were zero. The second scenario was Random Model. The sparse covariance graph was created by assigning to be non-zero with probability , independently of other elements.

Random Model Cliques Model
Method Root-Mean-Square Error Entropy Loss Root-Mean-Square Error Entropy Loss
Graphical Lasso
Table 1: Summary of simulation results for optimization problem (A-3). Sample size is chosen as . simulations were generated under Random Model and Cliques Modes. For each scenario, we present the empirical average and standard error of Root-Mean-Square Error and Entropy Loss.
Random Model Cliques Model
Method Execution Time Execution Time
Graphical Lasso
Table 2: Summary of simulation execution time for optimization problem (A-3). Sample size is chosen as . simulations were generated under Random Model and Cliques Modes. For each scenario, we present the empirical average and standard error of the execution time.

From the Tables 1 and 2, we can see that both Graphic-Lasso and Ledoit-Wolf optimization model had smaller mean square error and entropy loss compared to SPCOV in both scenarios. In addition, all the three algorithms were more accurate and robust in Random Model. In Figure 1, the different shades of color in the connectivity matrix indicated the correlation between each nodes and darker colors demonstrated stronger correlations. It can be seen that the estimated connectivity matrix obtained from Graphical Lasso and Ledoit-Wolf algorithms were in high accordance with the true covariance structure. In summary, simulations suggest that both Graphical Lasso and Ledoit-Wolf methods are competitive for identifying the sparsity structure of the simulated data.

The connectivity matrix of Random Model
The connectivity matrix of Cliques Model
Figure 1: Visualization of simulation results in Random Model

5 Real-data analysis

5.1 Dataset description

The EEG dataset was obtained from a multi-subject stroke experiment conducted at the University of California Irvine Neurorehabilitation Lab (PI: Cramer). During the experiment, participants sat in a chair facing a monitor in a single session. Their task was to make movements across different centers of each circle on the screen. To minimize the variability between individuals, the researchers measured the awake resting-state EEG for 3 min (EEG-Rest) at 1000 Hz prior to the motion task. Then, the measurement of each participant’s maximum arm movement speed was obtained, and a baseline assessment of motor skill task was recorded. During this procedure, EEG was measured (EEG-Test1). Later on, the participant was required to receive a practice block, followed by another test block. Finally, after three tests and two practice blocks were done, the EEG was obtained, which comprised of four scenarios – EEG-Test(1-3) and EEG-Rest. Figure 2 shows the general setup of this experiment.

Figure 2: Experiment setup

There are 16 subjects recorded in the entire dataset. As for each subject, each scenario of dataset consists of 160 trials, 1000 time points and 256 channels indicating different cortical regions. In this paper, subject “YUGR” is chosen. Figure 3 shows the EEG signal across trials in channel 6 from subject “YUGR”. The dataset consists of 160 epochs from EEG-Rest, 73, 74 and 63 trials from EEG-Test 1-3 respectively.

Figure 3: EEG signals across trials in channel 6

5.2 Modeling and results

To study the brain connectivity during the motion experiment, we applied the EEG data to the optimization problem (A-3). Motivated by the results from simulations, Graphical Lasso and Ledoit-Wolf algorithms were implemented to search for the solution. Figure 4 shows the connectivity matrix across different brain cortices from the two algorithms. It can be found that both of the matrices show a similar pattern in regrading to the correlation between cortices. In particular, the total cortices can be classified as four regions in which high association can be realized.

Figure 4: The connectivity matrix of the EEG data

6 Conclusion

We have developed a method to model brain connectivity through graphical models based on the framework of Copula Gaussian Model. To further address the optimization problem (A-3), Graphical Lasso, Ledoit-Wolf and SPCOV algorithms were implemented and compared. Simulations results show the benefit of using Graphical Lasso and Ledoit-Wolf in searching for the connectivity matrix. The proposed algorithms were also conducted on the real EEG data from the motion experiment. Results show the sparsity of the brain connectivity between cortices, which coincides with the results from previous literature.


  • Trudeau (1993) Trudeau, R.J. (1993). Introduction to Graph Theory. New York.
  • Bishop (2006) Bishop, C.M. (2006). Pattern Recognition and Machine Learning. Springer-Verlag.
  • Gao et al. (2016) Gao, X. & Shahbaba, B. & Fortin, N. & Ombao, H.(2016). Evolutionary State-Space Model and Its Application to Time-Frequency Analysis of Local Field Potentials. arXiv preprint arXiv:1610.07271.
  • Gao et al. (2018) Gao, X. & Shen, W. & Ombao, H. (2018). Regularized matrix data clustering and its application to image analysis. arXiv preprint arXiv:1808.01749.
  • Bullmore& Bassett (2011) Bullmore, E.T. & Bassett, D.S.(2011). Brain graphs: graphical models of the human brain connectome. Annual Review of Clinical Psychology. 7, 113–140.
  • Kemmer et al. (2015) Kemmer, P.B. & Guo, Y. & Wang, Y. & Pagnoni, G.(2015). Network-based characterization of brain functional connectivity in Zen practitioners. Frontiers in psychology. 6, 603.
  • Friedman et al. (2008) Friedman, J. & Hastie, T.& Tibshirani, R.(2008). Sparse inverse covariance estimation with the graphical lasso. Biostatistics. 9, 432–441.
  • Fiecas et al. (2010) Fiecas, M. & Ombao, H.& Linkletter, C.& Thompson, W.& Sanes, J.N.(2010). Functional Connectivity: Shrinkage Estimation and Randomization Test. NeuroImage. 40, 3005-3014.
  • Fiecas& Ombao (2011) Fiecas, M. & Ombao, H.(2011). The Generalized Shrinkage Estimator for the Analysis of Functional Connectivity of Brain Signals. Annals of Applied Statistics. 5, 1102-1125.
  • Dauwels et al. (2012A) Dauwels, J. & Yu, H.& Wang, X.& Vialatte, F.& Latchoumane, C.& Jeong, J.& Cichocki, A.(2012). Inferring Brain Networks through Graphical Models with Hidden Variables. Machine Learning and Interpretation in Neuroimaging. 7263, 194–201.
  • Achard& Bullmore (2007) Achard, S. & Bullmore, E.(2007). Efficiency and cost of economical brain functional networks. PLoS Computatinal Biology 3:e17.
  • Michelovannis et al. (2006) Michelovannis, S. & Pachou, E.& Stam, C.J.& Breakspear, M.& Bitsios, P.& et al.(2006). Small-world networks and disturbed functional connectivity in schizophrenia. Schizophrenia Research. 87, 60–66.
  • Stam (2004) Stam, CJ.(2004). Functional connectivity patterns of human magnetoencephalographic recordings: a “smallworld” network?. Neuroscience Letters. 355, 25–28.
  • Gong et al. (2008) Gong, G. & He, Y.& Concha, L.& Lebel, C.& Gross, DW.(2008). Mapping anatomical connectivity patterns of human cerebral cortex using in vivo diffusion tensor imaging tractography. C. Cerebral Cortex. 19, 524–536.
  • Dauwels et al. (2012B) Dauwels, J. & Yu, H.& Wang, X.(2012). Copula Gaussian Graphical Models with Hidden Variables. IEEE International Conference on Acoustics, Speech, and Signal Processing. Kyoto, Japan.
  • Bein& Tibshirani (2010) Bein, J. & Tibshirani, R.(2010). Sparse Estimation of a Covariance Matrix. Biomtrika.0, 1–24.
  • Rothman et al. (2008) Rothman, A & Levina, E. & Zhu, J.(2008). Sparse permutation invariant covariance estimation. Electronic Journal of Statistics.2, 494–515.
  • Butte et al. (2000) Butte, A & Tamayo, P. et al.(2010). Discovering functional relationships between RNA expression and chemotherapeutic susceptibility using relevance networks. Proceedings of the National Academy of Sciences of the United States of America.97, 12182–12186.
  • Rothman et al. (2009) Rothman, A & Levina, E. & Zhu, J.(2009). Generalized thresholding of large covariance matrice. Journal of the American Statistical Association.104, 177–186.
  • Rothman et al. (2010) Rothman, A & Levina, E. & Zhu, J.(2009). A new approach to Cholesky-based covariance regularization in high dimensions. Biometrika.97, 539.
  • Ledoit& Wolf (2012) Ledoit, O. & Wolf, M.(2012). Nonlinear shrinkage estimation of large-dimensional covariance matrices. The Annals of Statistics.40(2), 1024-1060.
  • Banerjee et al. (2007) Banerjee, O. & Ghaoui, L.E. & d’Aspremont, A(2007). Model selection through sparse maximum likelihood estimation. Journal of Machine Learning Research.101.
  • Dahlhaus (2000) Dahlhaus, R.(2000). Graphical interaction models for multivariate time series. Metrika.5, 157-172.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description