The human brain comprises approximately 1011 neurons, each of which makes about 103 synaptic connections. This huge number of connections between individual processing elements provides the fundamental scaffold for neuronal ensembles to become transiently synchronized or functionally connected . A similar complex network configuration and dynamics can also be found at the macroscopic scales of systems neuroscience and brain imaging . The emergence of dynamically coupled cell assemblies represents the neurophysiological substrate for cognitive function such as perception, learning, thinking . Understanding the complex network organization of the brain on the basis of neuroimaging data represents one of the most impervious challenges for systems neuroscience.
Several measures to evaluate at various scales (single cells, cortical columns, or brain areas) how the different parts of the brain communicate have been recently proposed. We can classify them, according to their symmetry, into two groups: symmetric and asymmetric measures. Symmetric measures, such as correlation, coherence, phase synchronization indexes (PLV, PLI), evaluate functional connectivity, i.e. statistically significant relationships between signals recorded from spatially remote neurophysiological events. On the other hand, the asymmetric ones are able to detect effective connectivity i.e. information flow (the influence one neuronal system exerts over another) , and therefore they reveal the direction of the interaction. Granger Causality (GC) belongs to this latter group.
2. Granger causality
2.1. Classical approach
A notion of causality, which is relevant for the experimental study of information transfer in neural systems, is commonly attributed to Wiener. He stated that, for two simultaneously measured signals, if one can predict the first signal better by incorporating the past information from the second signal than by using only information from the first one, then the second signal can be called causal to the first one . However, it was the Nobel Prize laureate Clive Granger who gave a mathematical formulation of this concept . He argued that, if X is influencing Y, then adding past values of the first variable (X) to the regression of the second one (Y) should improve the prediction of the latter.
To quantify such influence, the most straightforward approach consists in comparing the prediction of future values of one of the variables (say, X), by using two different models. The first one is a model where these values are forecast exclusively from the past history of this variable (e.g., a purely autoregressive (AR) model of a given order). The second one incorporates instead the past history of Y as well. Then, GC is defined by comparing the forecast errors of both models (see section 4.2), and Y is said to Granger-cause X if the error of the second model is significantly lower than that of the first one.
GC has the advantage of quantifying the information flow in the data, which is important to study the direction of the relationships between different parts of the brain. GC was first introduced in the field of economics 40 years ago. However, applications to neuroscience are rather more recent, see  for a review on applications to neurophysiology. One of the first studies using this concept investigated the existence of directional or causal interactions by analyzing local field potentials (LFPs) from the macaque inferotemporal cortex . This method was also applied to the LFP data recorded from two separate areas (primary and higher visual areas) of the cat visual cortex, to investigate the role of bottom-up and top-down interactions in a go/no-go task  or in a stimulus expectancy task . A frequency specific Granger causality measure was utilized in LFP recordings from somatosensory and motor cortices of macaque monkeys as they performed a motor maintenance in a visual discrimination task . In human, a time-variant Granger causality measure was applied to EEG signals from the standard color-word conflict Stroop task . A wide study CG’s advantages for neuroscience was done in .
3. Magnetoencephalography (MEG)
Ethical approval was granted by the local Ethics Committee. All data were analyzed anonymously. Subjects who underwent MEG recordings for research purposes had given written informed consent before participating.
Once subjects arrived at the laboratory, they got familiar with the MEG chamber in which the recording took place and were instructed about the experimental procedure. All the recordings took place during morning hours.
3.2. MEG system
MEG recordings were obtained using a 306-channel wholehead Elekta Neuromag® MEG system (Elekta Oy, Helsinki, Finland). The system has 102 magnetometers and 204 planar gradiometers in a helmet-shaped array covering the entire scalp. Subjects were seated inside a magnetically shielded room during MEG recordings (Vacuumschmelze GmbH, Hanau, Germany).
Eye movements were monitored by simultaneously recording an electrooculogram (EOG), with three Ag/Cl electrodes, two above and below the right eye and one at the right earlobe used as ground reference.
3.3. Artifact control
Participants were seated with their head in the MEG sensor helmet, which covers the entire head except the face. Four head position indicator coils (HPI) were placed on the scalp, appropriately spaced in the region covered by the MEG helmet. The locations of the nasion, two pre-auricular points, and the four HPI coils were digitized prior to each MEG study using a 3D-digitizer (FASTRACK; Polhemus, Colchester, VT) to define the subject-specific cartesian head coordinate system. 100-200 additional anatomical points were digitized on the head surface to provide a more accurate shape of the subject’s head. Once a subject was comfortably positioned in the MEG machine, short electrical signals were sent to the HPI coils to localize them with respect to the MEG sensor array. The data from the HPI coils were used to correct head movement during the session.
The data used in this study were acquired from 8 patients suffering from frontal focal epilepsy (FE), 8 patients suffering from generalized epilepsy (GE) and 8 healthy subjects (HS). Their ages range between (FE: 35.1 14, GE: 23.5 4, HS: 19 1).
3.5. Task and parameters
MEG data was acquired at a sampling rate of 1 kHz, with on-line band-pass filter of 0.10–330 Hz. Acquisition occurred in a single 20 min session of “resting state”; 10 min with eyes open looking to a cross on the screen followed by 10 min with eyes closed.
3.6. Removing artifacts
To correct the head position and the associated movement-related artifacts, a spatio-temporal signal space separation method (tSSS) with movement compensation (MC) was applied using the MaxFilter® software (Elekta Neuromag®, Elekta Oy, Helsinki, Finland). The use of tSSS is especially important for rejection of close-to-sensor artifacts.
After data reconstruction by signal space separation (referred as SSS), tSSS identifies artifacts by their correlated temporal behavior inside and outside the sensor helmet. The artifacts to be eliminated are thresholded by the quantitative level of this correlation determined by correlation limit of 0.95, and a correlation windows length of 10ms.
For each subject, artifact free epochs were then carefully selected by visual analysis Artifacts were typically due to (eye) movements, drowsiness or technical issues.
3.7. Used data
All the analysis were performed on 960 non-overlapping quasistationary segments (40 segments per subject) of 5000 ms free from eye or muscular artifacts, and epileptic activities (as, e.g., seizures or epileptic-like activity) and far (at least 20s) from recent epileptic discharges. The period of resting state with closed eyes was selected for the study. A downsampling to 500Hz was applied, thus obtaining segments of 2500 samples, as this has proven to be sufficient to detect clinically relevant differences in functional or effective connectivity in previous studies.
One important practical issue we would like to tackle here is that, in human MEG recording, not all types of magnetic sensors have the same sensitivity to deep brain sources. Magnetometers measure the overall magnitude of the magnetic field component approximately normal to the head surface; whereas planar gradiometers measure the difference of that field component at two adjacent locations. Describing MEG sensors in descending order of sensitivity to the depth of sources, magnetometers are most sensitive, followed by first-order axial gradiometers, second-order gradiometers and, finally, planar gradiometers . Hence, planar gradiometers have maximum sensitivity to sources directly under them, i.e., superficial cortical sources , which makes them less sensitive to artifacts and distant disturbances, and therefore are suitable for studying this case, as these kinds of epilepsy have a neocortical origin. However, for the analysis, data from both the planar gradiometers and the magnetometers, was used.
4. Data analysis
4.1. Data preprocessing
In agreement with previous findings, surrogate data tests revealed that less than 4% of interdependencies between the spontaneous brain activities were nonlinear . Thus, weighted brain networks were constructed by means of a definition of functional links based on Granger Causality (GC).
Application of GC requires that each time series is ‘covariance stationary’ (CS), i.e., that its mean and variance do not change over time. CS can be assessed in a rule-of-thumb way by examining the auto-correlation function. A non-CS time series will have an autocorrelation function that falls off slowly; a CS time series will have a sharply declining autocorrelation function. For this reason is that we take segments far enough from the spikes (at least 20s).
Prior to the causality analysis, we detrended the signals (subtracting the best-fitting line from each time series) and removed their temporal mean to provide a ‘zero-mean’ situation . Granger causality analysis was performed to each segment of data (Anil Seth Toolbox ). The order of the autoregressive (AR) model was set to 4 according to the Bayesian Information Criterion (BIC) for our particular data .
4.2. Granger causality
As commented above, for two simultaneously measured signals, if one can predict the first signal better by incorporating the past information from the second signal than using only information from the first one, then the second signal can be called causal to the first one .The simplest way of quantifying this is by using univariate vs. bivariate linear regression models of a signal x(t), to compare the forecast errors obtained incorporating (bivariate) or not (univariate) information from the past of y(t).
For the univariate autoregressive model (AR), we have:
where aij are the model parameters (coefficients usually estimated by least square method), P is the order of the AR model and ui, are the residuals associated to the model. Here, the prediction of each signal is performed only by its own past.
Accordingly, for the bivariate AR:
where the residuals uXY and uYX now depend on the past values of both signals, and their variance is:
where var(.) denote variance over time and is the prediction of x(t) by the past samples of values of x(t) and y(t).
Therefore, Granger causality from y(t) to x(t) (predicting X from Y) is:
GC ranges from 0 to infinity. The lower bound indicates that the incorporating information about the past of y(t) does not improve the prediction of x(t): . Accordingly,  is greater than 0, when the past of y(t) does improve the prediction of x(t): (Y G-causes X)
As indicated in the Introduction, GC belongs to the category of asymmetric connectivity indexes; therefore, it assesses effective rather than functional connectivity. However, as defined from AR models it is a linear parametric method, so it is only sensitive to linear correlations and depends on the order (P) of the autoregressive. (See  for an extension of the GC concept to assess nonlinear interdependencies)
4.2.1. Significance of GC values
If we calculated the GC index  for all possible pair of sensors of the same kind (e.g. magnetometers, planar gradiometers along one axis and planar gradiometers along the two orthogonal axis) we get three connectivity matrixes, where each element GCi->j measures the degree in which the signal in the sensor i Granger-causes that of the sensor j (i,j=1,..,102). However, it is necessary to determine whether the value of each given index is due to real causal dependence between the corresponding two sensors or rather it is not zero due to statistical fluctuations in the estimation of variance of the forecasting. This task can be easily accomplished because it is known that time-domain GC interaction is significant if the coefficients aij are jointly significantly different from zero. This can be established via an F-test on the null hypothesis that Aij are zero .
These have to be corrected for multiple comparisons, which we do here by applying the well-known ‘false discovery rate’ (FDR) method [21–23] at a desired significance threshold of 0.05. In all the cases where the (corrected) test failed to reject the null hypothesis, the corresponding GC index is set to zero.
Moreover, recent results show that correlations between magnetic fields sensors located at a distance less than 4 cm cannot distinguish between spontaneous activities of epileptic patients and control subjects . To reduce the influence of these spurious correlations between MEG signals, we have excluded the nearest sensors (separated less than 4 cm) from the computation of Granger causality values.
The final effective connectivity matrixes so obtained can be interpreted as the adjacency matrix of a directed weighted network. Hence, network analysis methods derived from graph theory  could be applied to assess the most significant features of the brain connectivity networks in the groups of subjects described before (i.e. GE, FE and HS), during the resting state task with closed eyes [26,27]. And from that, study their differences.
4.3. Network parameters
To characterize the network structure of brain activity of each group of subjects, we evaluated a list of measures for directed weighted graphs . In this approach, MEG sensors were considered as vertices (nodes) and the GC values between sensors as edge weights (links). The edge weight represents the strength of the connection between the vertices.
Here some notation is explained (see  for details):
N is the set of all nodes in the network, and n is the number of nodes.
L is the set of all links in the network, and l is number of links.
(i, j) is a link between nodes i and j, where i,j ∈ N.
Links (i, j) are associated with connection weights wij.
aij is the connection status between i and j: aij = 1 when link (i, j) exists (when i and j are neighbors); aij = 0 otherwise (aii = 0 for all i).
We compute the number of links as l =∑i,j∈N aij (to avoid ambiguity with directed links we count each undirected link twice, as aij and as aji).
The sum of all weights in the network is lw, and it is computed as lw =∑i,j∈N wij. Henceforth, we assume that weights are normalized, such that 0 ≤ wij ≤ 1 for all i and j.
We focused on the following global parameters, to be explained henceforth: the average degree (section 4.3.1), strength (section 4.3.2) and two measures of segregation, namely the clustering coefficient (section 4.3.3) and modularity (section 4.3.4).
The degree (D) of a node is the number of links connected to it. In directed weighted networks, we distinguish between the ‘in degree’, which is the number of links that arrive to the node, and the ‘out degree’, which is the number of links that go out from the node.
(Directed) in-degree of i, kini =∑j∈Naji
(Directed) out-degree of i, kouti =∑j∈Naij
The global degree of a network is the average of all its nodes’ degree. The mean network degree is most commonly used as a measure of density, or the total “wiring cost” of the network.
The weighted variant of the degree, sometimes termed the strength, is defined as the sum of all neighboring link weights.
4.3.3. Clustering coefficient
The clustering coefficient (C) describes the likelihood that neighbors of a vertex are also connected. It quantifies the tendency of network elements to form local clusters. We used the directed and weighted equivalent of this measure to characterize local clustering :
Here and in the following equation, the arrow next to the indexes denotes that they are calculated for directed networks. However, for the sake of clarity in the text, we will omit this arrow in the name of the indexes henceforth.
Modularity (Q) quantifies how a network can be optimally divided in subgroups or modules. We used a modification for directed weighted networks by 
Roughly speaking, the greater the value of Q, the more modular a network is, i.e., the greater the density of within-cluster connections as compared to the between-cluster ones.
4.4. Statistical tests
4.4.1. Between groups comparisons
We performed a Kruskal Wallis test to compare network parameters from the three groups (FE, GE, HS). This test compares the medians of the samples in each group, and returns the p value for the rejection of the null hypothesis that all samples are drawn from the same population (or equivalently, from different populations with the same distribution). The Kruskal-Wallis test is a nonparametric version of the classical one-way ANOVA, and an extension of the Wilcoxon rank sum test to more than two groups.
If the p value is close zero, this suggests that at least one sample median is significantly different from the others. Here, we took p<0.05 as the critical value to determine whether the null hypothesis can be rejected.
4.4.2. Post-hoc between (two) groups comparisons
In those cases where the Kruskall-Wallis to check was significant, we further analyzed pairwise difference between any two groups by means of a two-sided rank sum test (Wilcoxon test), test the null hypothesis that the network parameters obtained for each group of subjects are independent samples from identical continuous distributions with equal medians, against the alternative that they do not have equal medians. This test is equivalent to a Mann-Whitney U-test. Differences were again considered significant if p<0.05.
Differences in the global network parameters between the three groups (FE, GE and HS):
5.1. Differences among the planar gradiometer networks
The interdependence matrixes obtained from the planar gradiometers were sparse, i.e., they presented only a few significant GC values (effective links) between the corresponding neocortical sources. Differences among groups for the planar gradiometers were found in modularity (see figure 1).
Patients suffering from frontal focal epilepsy, present a higher value of Q as compared to those suffering from generalized epilepsy and healthy subjects. Differences in magnetometers
5.2. Between-group differences in brain networks recorded from magnetometers
As can be seen from figure 2, the FE group presents a lower value of C than those from the GE and the HS group.
Also, for the average degree D (figure 3) patients from the FE group, present a higher value of global network degree as compared to GE patients and healthy controls.
Moreover, for the average strength S (figure 4), FE patients showed the highest values of the index as compared to GE patients and healthy controls.
Finally, for the number of clusters (figure 5) healthy subjects present a lower number of clusters as compared to both group of epileptic.
The assessment of functional brain networks from multivariate fMRI, EEG or MEG data has become a very popular line of research nowadays [25,31]. This is most likely because the parameters derived from this approach are relatively easy to calculate and interpret in neurophysiological terms. In fact, one only needs, on the one hand, a measure of functional or effective connectivity between any two signals and, on the other, a way of estimating the significance of such a measure. Once you are equipped with these two tools, the characterization of the corresponding brain network is only a few calculations ahead.
Yet even these apparently straightforward steps (e.g., estimation of bivariate connectivity index and their significance) are deceptively simple, and their careless, mechanical application may result in misleading results and/or false conclusions. Common sources of error are the misuse of a given connectivity index, failing to assess its significance in a proper way or not taking into account what the data we have recorded are actually measuring. This latter issue is specially significant when dealing with extra-cranially recorded signals (whether EEG or MEG data), in which it is well-known that signals in different sensors are often measuring the activity of the same deep brain source, this resulting in a statistically significant relationship between these signals that has very few (if any) to do with the existence of connectivity between the sensors.
This important question can be tackled in different ways. One obvious attempt consists in trying to reconstruct, from the recorded data, the neurological sources of activity before applying the network connectivity approach. But the so-called inverse problem of source reconstruction is known to be ill-posed, in the sense that it is underdetermined, as the number of sources is much larger than that of the measures, and has not a single solution. Another possibility includes the use of carefully devised indexes , which are robust against common source effects but which, unfortunately, are not without their own particular problems . In this work, we have dealt with the question in a different way, by taking advantage of the complementary information provided by two kinds of sensors (magnetometers and planar gradiometers) in a modern MEG recording device. The key idea consists in the fact that whereas magnetometers are sensitive to all sources (whether cortical or deep ones) in the direction normal to the head surface where they record, gradiometers are only sensitive to cortical sources right behind the recorded surface , and therefore unlikely to be affected by the common source problem. Thus, comparing the results from the magnetometers and the planar gradiometers recorded simultaneously from the same subject, we may get insight into the possible effect that common sources has on the functional brain network assessed by the former sensors, which are by far the most common sensors used in MEG literature.
Here, we have applied this approach to study brain networks of effective connectivity from planar gradiometers and magnetometers recorded during rest with closed eyes in three groups of subjects, one of healthy controls and two of patients suffering from different kinds of epilepsy: generalized and frontal focal. The connectivity index used was the GC index, whose significant was assessed by a careful combination of F-test estimation and FDR statistical test to deal with the multiple comparison problem. And we think we might rightly conclude that the results obtained justify our use of what we may call this multi-sensor approach. Indeed, thanks to it, two important questions have come out. First, if we were to analyze the data from the magnetometers alone, we would have come to the conclusions that the functional networks from the three groups are different in a number of parameters, whereas these networks, when only the cortical sources are recorded, do not present differences among the three groups in any of these parameters. Second, networks of cortical sources do present differences in modularity between healthy and epileptic groups, as well as between the two epileptic patients.
Thus, by analyzing the information provided by the combined use of both sensors, we have a greater insight into the brain networks of the three groups (and the differences among them). The information from the planar gradiometers suggests that connectivity between neocortical sources is sparse, and the corresponding networks are only weakly modular (low Q), but the intra-cluster connectivity of the epileptic groups are greater than that of healthy subjects, with the FE group in turn being significantly more modular than the GE group. This result is very interesting, as it is well-known that the greater the intra-cluster connectivity, the easier it is for a network to become fully synchronized . Thus, it would indicate that in both types of epilepsy, the network of neocortical sources is more prone to become (pathologically) synchronized that in normal subjects.
When we review the information from the magnetometers, we can get information about deeper brain sources. The first interesting result is that the networks constructed from magnetometers are more densely connected that those from the gradiometers. Besides, these networks are scarcely modular in any of the groups (as the number of clusters in Figure 5 is close to the number of sensors) and present a higher number of (stronger) links (higher C, D and S as shown in figures 2, 3 and 4, respectively) in the epileptic groups as compared to healthy subjects. Unfortunately, it is difficult to interpret these results simply in terms of the corresponding connectivity between deep sources, as the possibility of the same source being recorded in different magnetometers cannot be ruled out. Yet rather than whether these results are a consequence of more active deep sources or more interconnected ones, what is really important is the fact that network theory as applied to MEG data is able to disclose changes in neurological activity of epileptic subjects as compared to healthy subjects even when no epileptogenic activity is apparent in the raw data.
6.1. Future work
The results presented above are certainly interesting, not only because of their methodological and neurological implications, but also because they pave the way for future studies in the same line. One obvious question to further investigate from these data is whether the between-groups differences found in various global parameters in both planar gradiometers and magnetometers can be explored in terms of their (possibly local) origin. Namely, whether they are due to local differences between the corresponding networks, which are detected by these global parameters but can be also characterized topologically. Indeed, it is tempting to hypothesize that the FE group, whose epilepsy has a focal origin, presents differences with both the GE and the HS group precisely in the frontal sensors. Another important question that we should elucidate is whether differences are present at all the frequencies of the signal or are rather constrained to a certain frequency band (as found in former works ), but detected here by the GC index, which operates on the broadband signal in the time domain. One last issue of practical interest that we are currently studying is the possible usefulness of the differences in network connectivity patterns to assist in the classification (diagnosis) of recorded subjects using machine learning algorithms. This latter application is potentially very important from the clinical point of view. Additionally, it would also allow circumventing the problem of determining if the patterns of GC index (or any other connectivity index for that matter) are due to true connectivity or are the result of volume conduction of deep brain sources, which are reflected in many sensors at the same time. As long as these patterns are different in each group, they would be useful to detect deviations from healthy condition. Here also, the fact that the information from the two type of magnetic sensors available (magnetometers and planar gradiometers) is not redundant but complementary, speaks clearly in favor of analyzing both sensors at the same time, whenever possible.
We hope that this line of research will continue to provide further insight into the pattern of connectivity networks in health and disease and its possible application as an additional tool for diagnostic purposes in different neurologic pathologies.
Guiomar Niso, Ernesto Pereda and Francisco del Pozo acknowledge the financial support of the Spanish Ministry of Economy and Competitiveness through the grants TEC2012-38453-CO4-01 and -03. Fernando Maestú, María Gudín and Sira Carrasco acknowledge the financial support of the Castilla La Mancha Sanitary System through grant PS09/00450. Guiomar Niso has also received the financial support of the Spanish Ministry of Education and Science through the FPU grant AP2008-02383.