International Journal of Bioelectromagnetism Vol. 8, No. 1, pp. 2-19, 2006.

www.ijbem.org

Estimate of Causality Between Cortical
Spatial Patterns During Voluntary Movements
in Normal Subjects

Hovagim Bakardjian^a, Andrzej Cichocki^a, Febo Cincotti^b, Donatella Mattia^b, Fabio Babiloni^c,b, Maria Grazia Marciani^b, Fabrizio De Vico Fallani^d,b, Fumikazu Miwakeichi^e, Yoko Yamaguchi^f, Pablo Martinez^a, Serenella Salinari, Andrea Tocci^g,and Laura Astolfi^b,g

^a Laboratory for Advanced Brain Signal Processing Riken, Brain Science Institute, Japan
^b IRCCS "Fondazione Santa Lucia", Rome, Italy
^c Dipartimento di Fisiologia umana e Farmacologia, Rome, Italy
^d Centro di Ricerca de "La Sapienza" per l'Analisi dei Modelli e dell'Informazione nei Sistemi Biomedici, Rome,Italy
^f Laboratory for Dynamics of Emergent Intelligence, Riken, Brain Science Institute, Japan
^e Department of Medical System Engineering, Faculty of Engineering, Chiba University, Japan
^g Dept. of Informatica e Sistemistica, Univ. of Rome “La Sapienza”, Rome, Italy

Correspondence: Laura Astolfi, IRCCS “Fondazione Santa Lucia”, Rome, Italy,
via Ardeatina 306, 00100, Rome, Italy 
E-mail: laura.astolfi@uniroma1.it phone +39 06 51501466

Abstract. Different non-invasive brain imaging techniques are presently available to provide images of the human functional cortical activity, based on hemodynamic metabolic or electromagnetic measurements. However, static images of brain regions activated during particular tasks do not convey the information of how these regions communicate with each other. Cortical connectivity estimation aims at describing these interactions as connectivity patterns which hold the direction and strength of the information flow between cortical areas. So far, the causality between brain signals has been assessed by using the time-varying information derived from hemodynamic or electromagnetic signals recorded at the scalp or cortical level. However, the causality estimation from these brain functional waveforms will depict a single pattern of connectivity involving several brain areas, for each time segment or in every frequency band analyzed. Since it is well known that the brain does not produce any “single waveform” but rather engages several distributed cortical areas in order to process information, a question arose about the appropriateness of the estimation of the functional connectivity between waveforms. In particular, the question is whether instead of estimating the causality between single waveforms derived from the different cortical or scalp areas it is possible to estimate the causality between “spatial patterns of brain cortical activations”. In fact, it is reasonable to pose the question if it could be more interesting to estimate the causality (in the sense of the Granger) between the activation of distributed cortical systems or just observe the causality between isolated waveforms. In this report, we attempted to estimate the causality between distributed cortical systems during the execution of simple movements in a group of normal healthy subjects. To estimate the causality between the spatially distributed patterns of cortical activity in the frequency domains we applied a series of processing steps on the recorded EEG data. From the high-resolution EEG recordings we estimated the cortical waveforms for the regions of interest (ROIs), each representing a selected sensor group population. The solutions of the linear inverse problem returned a series of cortical waveforms for each ROI considered and for each trial analyzed. In each subject, the cortical waveforms were then subjected to an Independent Component Analysis (ICA) pre-processing. The independent components obtained by the application of the ThinICA algorithm (which combined second- and fourth-order statistics to achieve spatial and temporal

decorrelation and mutual independence of extracted latent (hidden) components) were then further processed by the Partial Directed Coherence algorithm, in order to extract the causality between the spatial cortical patterns of the estimated data. Such couples of cortical patterns were obtained for each one of the three frequency bands employed here; theta, alpha, and beta. The estimated cortical patterns indicated the involvement of a large network of parietal and premotor areas in the beta band, as well as the existence of a distributed network in the theta and alpha frequency bands that was driving consistently the right primary motor and supplementary motor areas. These results are the first to demonstrate the involvement of a group of cortical areas that “causes” the activation of a second group of cortical areas for a simple motor task.

Keywords:ThinICA; Distributed current density estimates; Brodmann areas; inverse problem; high-resolution EEG, Functional connectivity, Partial Directed Coherence

1. Introduction

Nowadays, different non-invasive brain imaging techniques are presently available to provide images of the human functional cortical activity, based on hemodynamic (functional Magnetic Resonance Imaging, fMRI), metabolic (Positron Emission Tomography, PET) or electromagnetic (Electroencephalography, EEG and Magnetoencephalography, MEG) measurements. However, static images of brain regions activated during particular tasks do not convey the information of how these regions communicate to each other. In fact, the concept of brain connectivity is viewed as central for the understanding of the organized behavior of cortical regions beyond the simple mapping of their activity [Lee et al., 2003, Horwitz, 2003, David et al., 2004]. This organization is thought to be based on the interaction between different and differently specialized cortical sites. Cortical connectivity estimation aims at describing these interactions as connectivity patterns which reflect the direction and strength of the information flow between cortical areas. To achieve this, several methods have been already applied on data gathered using both hemodynamic and electromagnetic techniques [Buchel and Frison, 1997, Gevins et al., 1989; Urbano et al., 1998, Brovelli et al., 2002; Taniguchi et al. 2000]. Two main definitions of brain connectivity have been proposed over these years: functional and effective connectivity [Friston, 1994]. While functional connectivity is defined as the temporal correlation between spatially remote neurophysiologic events, the effective connectivity is defined as the simplest brain circuit which would produce the same temporal relationship as observed experimentally between cortical sites. As for the functional connectivity, the several computational methods proposed to estimate how different brain areas are working together typically involve the estimation of some covariance properties between the different time series measured from the different spatial sites during motor and cognitive tasks studied by EEG and fMRI techniques [Gevins et al., 1989; Urbano et al., 1998; Gerloff et al., 1998; Jancke et al., 2000]. So far, the estimation of functional connectivity on EEG signals has been addressed by applying either linear and non-linear methods which can both disclose the direct flow of information between scalp electrodes in the time domain, although with different computational demands [Nunez, 1995; Clifford et al., 1987; Inouye et al., 1995; Stam et al., 2002; Stam et al., 2003, Tononi et al., 1994; Quian Quiroga et al., 2002]. In addition, due to the evidence that important information in the EEG signals are often coded in the frequency rather than time domain (reviewed in [Pfurtscheller and Lopes da Silva, 1999]) attention has been focused on detecting frequency-specific interactions in EEG or MEG signals by analyzing the coherence between the activity of pairs of structures (Mizuhara et al., 2005). Coherence analysis does not have, however, a directional nature (i.e. it just examines whether a link exists between two neural structures, by describing instances when they are in synchronous activity), and it does not provide directly the direction of the information flow. In this respect, a multivariate spectral technique called Directed Transfer Function (DTF) was proposed [Kaminski et al., 2001] to determine the directional influences between any given pair of channels in a multivariate data set. This estimator is able to characterize at the same time the direction and the spectral properties of the brain signals, requiring only one multivariate autoregressive (MVAR) model to be estimated from all the EEG channel recordings. The DTF technique has been recently demonstrated to rely on the key concept of Granger causality between time series [Granger, 1969], according to which an observed time series x(n) causes another series y(n) if the knowledge of x(n)’s past significantly improves prediction of y(n); this relation between time series is not reciprocal, i.e. x(n) may cause y(n) without y(n) necessarily causing  x(n). This lack of reciprocity allows the evaluation of the direction of information flow between structures.

So far, the causality between brain signals has been assessed by using time varying information derived from hemodynamic or electromagnetic signals recorded at the scalp level [Kaminski et al., 2001] or estimated at the cortical level [Babiloni et al., 2005; Astolfi et al., 2005]. However, the causality estimation from these brain functional waveforms will depict a single pattern of connectivity involving several brain areas, for each time segment or in every frequency band analyzed. Since it is well known that the brain does not produce any “single waveform” but rather engages several distributed cortical areas in order to process information, a question arose about the appropriateness of the estimation of the functional connectivity between waveforms. In particular, the question is whether instead of estimating the causality between single waveforms derived from the different cortical or scalp areas, it is possible to

estimate the causality between “spatial patterns of cortical brainactivations”. In fact, it is reasonable to pose the question if it could be more interesting to estimate the causality (in the sense of the Granger) between the activation of distributed cortical systems or just observe the causality between isolated waveforms. In this report, we attempted to estimate the causality between distributed cortical systems during the execution and the imagination of simple movements in a group of normal healthy subjects. To estimate the causality between the spatial distributed patterns of cortical activity in the frequency domains we applied a series of processing steps on the recorded EEG data.

First, we estimated the cortical activity from the EEG recordings by using realistic head and cortical models for each subjects. This was obtained by considering the cortical activity occurring in the Brodmann areas used to segment the cortical models used. The cortical activities were estimated by using the solutions of the EEG linear inverse problems as described previously [Babiloni et al., 2005; Astolfi et al., 2005].

Furthermore, we applied Independent Component Analysis (ICA) pre-processing on the cortical waveforms derived from the Brodmann areas for each subject. The application of the to the cortical waveforms derived from the Brodmann areas returned a series of basic spatial patterns of activations as well as the temporal variation of these patterns along the estimated cortical waveforms.

The key point of all this processing is the estimation of causality between the temporal independent components calculated for the cortical data by using the DTF or Partial Directed Coherence (PDC, Baccalà and Sameshima, 2001). The application of the DTF or PDC between two different independent components estimated from the computed cortical waveforms returned an estimation of the causality between the two distributed cortical activation patterns. Each cortical activation pattern represented a spatial component and corresponded to an analysed independent component. Hence, the causality of the activation of the distributed cortical areas occurred not at the level of single waveforms, but rather at the level of a coordinated series of cortical areas, as described by the spatial independent components obtained by the ICA.

In this study, we propose to estimate the causality of the cortical connectivity patterns by exploiting the combined use of ICA and DTF/PDC techniques applied to high-resolution EEG signals which exhibit a higher spatial resolution than conventional cerebral electromagnetic measurements made over or outside of the scalp. The high-resolution EEG technique includes the use of a large number of scalp electrodes, realistic models of the head derived from structural magnetic resonance images (MRIs), and advanced processing methodologies related to the solution of the linear inverse problem. These methodologies allow the estimation of cortical current density from sensor measurements [Babiloni et al., 2000; Grave de Peralta and Gonzalez Andino, 1999; Pascual-Marqui, 1995; He et al., 2002]. Subsequently, a novel combination of ICA and DTF/PDC methods was applied to the cortical estimates obtained from high-resolution EEG data related to a simple movement task in humans.

2. Material and Methods

2.1. The estimation of the cortical activity from high-resolution EEG recordings

High-resolution EEG recordings

Six right-handed healthy subjects participated in the study. Informed consent was obtained from each subject after explanation of the study, which was approved by the local institutional ethics committee. For the EEG data acquisition, subjects were comfortably seated on a reclining chair, in an electrically shielded, dimly lit room. All the subjects gave their informed consent to participate in the study. In the present event-related experimental design, we adopted a simple motor task consisting of repetitive self-generated overt executions (control subjects) of the right foot dorsal flexion at the ankle. The absence of external cues was chosen in order to avoid that any part of the observed EEG task-induced activity could be related to the perception or processing of pacing stimuli per se. A 58-channel EEG system (BrainAmp, Brainproducts GmbH, Germany) was used to record electrical potentials by means of an electrode cap with sensors placed according to the extended 10-20 international system. Structural MRIs of the subject’s head was taken with a Siemens 1.5T Vision Magnetom MR system (Germany). The EEG was sampled at 200 Hz, and 100 trials were recorded for each subject. Figure 1 presents the experimental subjects for the normal group recorded with the high-resolution EEG cap.

Applying the tools for the estimation of cortical activity and connectivity

We estimated the cortical activity from the high-resolution EEG recordings, by using realistic head models and a cortical surface model with an average of 5,000 dipoles, which were uniformly distributed. The estimation was obtained by the application of a linear inverse procedure [Rolando Grave de Peralta and Gonzalez Andino, 1999; Babiloni et al., 2005]. Cortical activity was then estimated in ROIs generated by the segmentation of the Brodmann areas (B.A.) on the accurate cortical model used. Bilateral ROIs considered in this analysis were: the primary motor areas for the foot (MIF) and the lip movement (A4_Lip), the proper supplementary motor area (SMAp),

the standard premotor area (A6), the posterior parietal areas (A7) and the cingulated motor area (CMA). The labels of the cortical areas also contain a postfix characterizing the considered hemisphere (R, right, L, left). Such ROIs were segmented on the basis of Talairach coordinates and anatomical landmarks available. Figure 2 presents the cortical areas as obtained for the realistic head models generated for each subject. It is possible to note the different positions of the same cortical areas on the cerebral surface between the subjects.

Figure 1. The normal subjects recorded in this experimental study, with the high-resolution EEG cap.

Figure 2. The cortical regions of interest (ROIs) employed in this study for the normal population investigated. Each ROI is represented with a different color, and the used color scheme is common across the different subjects. Note that the Cyngulate Motor Areas, located in the mesial central part of the cortical surface, are hidden in the interhemispheric scissure.

For each time point of the recorded EEG we solved the linear inverse problem, estimated the magnitude for each one of the thousands of dipoles used for cortical modelling. Then, we computed the average of the magnitude of such dipoles in each ROI considered, for each time point considered. The resulting cortical waveforms, one for each predefined ROI, were then subjected to the ICA analysis to reveal their independent components. Since the cortical waveforms were obtained for each single trial and recorded for each experimental subject, they were concatenated to each other, after detrending, for consecutive analysis as described in details below.

2.2 The estimation of the functional connectivity by DTF and PDC

Multivariate Methods for the Estimation of Connectivity

Let Y(t) be a set of cortical waveforms, obtained from several cortical regions of interest (ROI) as described in detail in the following way:

where t refers to time and N is the number of cortical areas considered.

Let’s assume that the following MVAR process is an adequate description of the data set Y:

where Y(t) is the data vector in time, E(t)=[e₁(t), …, e_N]^T is a vector of multivariate zero-mean uncorrelated white noise processes, Λ(1), Λ(2), … Λ(p) are the NxN matrices of model coefficients and p is the model order. In the present study, p was chosen by means of the Akaike Information Criteria (AIC) for MVAR processes [Akaike 1974] and was used for MVAR model fitting to simulations, as well as to experimental signals. It has been noted that, although the sensitivity of MVAR performance depends on the model order, small model order changes do not influence results (Franaszczuk et al, 1985; Babiloni et al, 2005).

The Directed Transfer Function

The Directed Transfer Function, representing the causal influence of the cortical waveform estimated in the j-th ROI on that estimated in the i-th ROI is defined (Kaminski and Blinowska, 1991] in terms of elements of the transfer matrix H:

In order to compare the results obtained for cortical waveforms with different power spectra, a normalization can be performed by dividing each estimated DTF by the squared sums of all elements of the relevant row, thus obtaining the so-called normalized DTF [Kaminski and Blinowska, 1991]:

where γ_ij(f) expresses the ratio of influence of the cortical waveform estimated in the j-th ROI on the cortical waveform estimated in the i-th ROI, with respect to the influence of all the estimated cortical waveforms. Normalized DTF values are in the interval [0, 1], and the normalization condition

is applied.

From the transfer matrix, we can calculate power spectra S(f). If we denote by V the variance matrix of the noise E(f), the power spectrum is defined by

where the superscript H denotes transposition and complex conjugate.

From S(f), standard coherence can be computed as:

Coherence measures express the degree of synchrony (simultaneous activation) between areas i and j.

Partial Directed Coherence

Partial coherence is another estimator of the relationship between a pair of signals, describing the interaction between areas i and j when the influence due to all N-2 time series is discounted. It is defined by the formula:

where M_ij(f) is the minor obtained by removing i-th row and j-th column from the spectral matrix S.

In 2001, Baccalà proposed the following factorization:

(12)

where Λ _n ( f ) is the n-th column of the matrix Λ ( f ). This led to the definition of Partial Directed Coherence (PDC) [Baccalà and Sameshima, 2001]:

(13)

The PDC from j to i, π_ij(f), describes the directional flow of information from the activity in the ROI s_j(n) to the activity in s_i(n), whereupon common effects produced by other ROIs s_k(n) on the latter are subtracted, leaving only a description that is specifically from s_j(n) to s_i(n).

PDC values are in the interval [0, 1], and the normalization condition

is verified. According to this condition, π_ij ( f ) represents the fraction of the time evolution of ROI j directed to ROI i, compared to all of j’s interactions with other ROIs.

For both DTF and PDC, high values in a frequency band represent the existence of an influence between any given pair of areas in the data set. However, an important difference is that PDC does not involve the inversion of matrix Λ. This leads to several points. In fact, an analysis of the definition of DTF reveals that, due to this matrix inversion, it is a linear combination of both the direct influence from one area to the other and the influence mediated by other areas, along various cascade pathways [Kaminski et al., 2001]. This becomes immediately clear from an example: given a three-region model, the non-normalized DTF from area 1 to area 2 is:/

θ212 ( f ) = | H 21 ( f )| 2 =   | Λ 21 ( f ) Λ 33 ( f ) – Λ 31 ( f ) Λ 23 ( f ) | | Λ ( f ) | 2

(15)

From this formula it can be noted that even if the direct influence from area 1 to area 2, Λ₂₁( f ), is zero, θ₂₁( f ) may still be different from zero, since there is an influence from 1 to 3 (Λ₃₁( f ) and from 3 to 2 (Λ₂₃( f ) ). The link between 1 and 2 will be indicated by DTF as a causal pathway if all the causal influences along the way are non-zero.

PDC, due to the lack of the matrix inversion, behaves differently. It indicates only the existence of a direct causal influence from area 1 to area 2. If no direct influence exists, PDC₂₁ is virtually zero.

2.3 The estimation of the causality between the distributed cortical patterns

As described above, the estimation of the cortical activity from the EEG recordings returns a cortical current density waveform for each ROI considered, in each trial analysed. In order to apply the ICA to the cortical current density waveforms, the following steps have been applied:

The procedure, as detailed above, will return an Anorm matrix that describes the mixing of the independent spatial cortical patterns for the gathered data, as well as a matrix Yscal revealing their temporal behavior. Hence, each i-th component of the Yscal matrix is related to the time-varying presence of the i-th spatial pattern described by the i-th column of the Anorm matrix in the gathered EEG data, represented by the X matrix.

Now, in order to assess the possible causality relations between the cortical patterns, we estimated the connectivity values between the different independent components described by the matrix Yscal. The connectivity values were obtained by using the DTF or PDC algorithm. In the following, the results will be presented as obtained by the PDC algorithm, while the results from the application of the DTF algorithm were similar to those described here.

2.4 The estimation of the causality between the patterns of the cortical activity

As described above, the application of the PDC algorithm to the Yscal matrix returned a series of causality patterns between the different independent components, each one related to a particular spatial pattern. The causality patterns between the independent components were considered only if statistically significant, in agreement with the procedure already described previously [Astolfi et al., 2005]. A statistical connection between the i-th and the j-th components of the Yscal matrix (represented as Yi –> Yj) means that the series of cortical ROIs involved in the i-th spatial pattern of the Anorm matrix will cause an activity in the series of cortical ROIs involved in the j-th spatial pattern of the same matrix Anorm. In the sense of the Granger theory, the inclusion of the Yi independent components (i.e. the cortical areas of the i-th independent component) improves the predictions of the time series of the Yj independent components in the multivariate autoregressive model.

For each subject analyzed, and for each one of the four frequency bands investigated (theta (4-7 Hz), alpha (8-12 Hz), beta (12-30 Hz) and gamma (above 30 Hz)), only the four highest connectivity links were used for the successive analysis. Such connectivity links were those with the highest statistical significance with respect to the random values of the DTF or PDC computed by using a shuffling procedure [Kaminski et al., 2001; Astolfi et al., 2005; Babiloni et al., 2005].

Hence, for each frequency band and every subject, we obtained a series of four cortical spatial patterns that “caused” or drove other four cortical spatial patterns during the execution of the task. In particular, these patterns were related to the preparation of the task for a time period of 1.5 seconds.

Since the independent components are not ordered between subjects, a possible problem arose when a comparison of these spatial patterns between the subjects have to be performed, in order to extract inferences related to the group behavior. In fact, it is well known that the numbering of the spatial components is not consistent between different subjects, i.e. the component number 4 for the subject k-th may not be the same as the component number 4 for the subject j-th, and so forth. Then, it is interesting to have a tool able to couple the independent components between subjects on the basis of their spatial patterns. This is important since the goal is to build a set of couples of cortical spatial patterns (one that “drives” and the other that it was driven) that are common for all the investigated populations. In order to obtain such “average” cortical pattern a series of operations have to be performed. The approach pursued in this study can be described in the following way:

Following this procedure, it is possible to obtain four “average” couples of spatial patterns for each frequency band. Such “average” cortical patterns correspond to the time period of the preparation for the actual movement in the normal healthy subjects.

2.5 Independent Component Analysis (ICA) and the ThinICA algorithm

ICA is a process which can extracts a new set of statistically independent components represented by the n-dimensional vector y(t) = W_x(t) from exploratory (observed) input data represented by the m-dimensional vector x(t) (t = 1, 2,… , N). The extracted components correspond to estimates of hidden or latent variables in the data sometimes called sources.  This process assumes that a time series  x(t) has an embedded mixing process of the form x(t) =A s(t), where A denotes an unknown mixing matrix and s(t) is a vector representing unknown hidden (latent) variables or sources. ICA can be considered as a demixing or a decomposition process which is able to recover the original sources, i.e., y(t) = ŝ(t) through the linear transformation y(t) = W_x(t). The fact that two random variables are uncorrelated does not also imply that they are independent. This fact is lost in using other methods such as Principal Component Analysis (PCA). The ICA approach seeks to find such independent directions through maximization of a suitable cost function called sometimes contrast function which is a measure of statistical independence. Such functions can be maximized or minimized using various optimization methods, including artificial neural networks.

Independent component analysis can be considered as an extension of principal component analysis (PCA). In PCA, the input data x(t) is decorrelated to find the components that are maximally uncorrelated according to second-order statistics. PCA gives orthogonalized and normalized outputs according to the second-order statistics by minimizing the second-order moments. The principal components can still be dependent however. The problem of independent component analysis or blind source separation of sources mixed instantaneously can be defined as follows. Let’s assume that we have available to us a set of multivariate time series {x_t(t)} (i = 1, 2, ... , m). We also assume that these time series, for example corresponding to individual EEG electrodes, are the result of an unknown mixing process defined by

or equivalently in compact matrix form X(t) = A S(t), (t = 1, 2, ... N), where A is an unknown mixing matrix sized m by n, and s(t) = [s₁(t), s₂(t), ... , s_n(t)]^T are hidden (latent) components called the sources. We seek to estimate the unknown sources s_j(t) using only the observed data vector x(t) = [x₁(t), x₂(t), ..., x_m(t)]^T. The problem is to find a demixing or separating matrix W such that y(t) = W x(t) estimates the hidden independent components. It is possible that there could be a different numbers of sensors than sources, that is, A may not be square. If it is assumed that the number of sources (hidden components) is the same as the number of time series or observed inputs n, then A is a square (n by n) matrix. If W = A^-1, then y(t) = s(t), and perfect separation occurs.

In practice, the optimal y will be some permutated and scaled version of s, since it is only possible to find W such that WA = PD where P is a permutation matrix and D is a diagonal scaling matrix. In general, the ICA of a random vector x(t) is obtained by finding a n by m, (with m ≥ n), full rank separating (transformation) matrix W such that the output signal vector y(t) = [y₁(t), y₂(t), ..., y_n(t)]^T, (independent components) estimated by y(t)=W x(t), are as independent as possible.

Compared with the principal component analysis (PCA), which removes second-order correlations from observed signals, ICA further removes higher-order dependencies. Statistical independence of random variables is a more general concept than decorrelation. Overall, we can state that random variables y_i(t) and y_j(t) are statistically independent if knowledge of the values of y_i(t) provides no information about the values of y_j(t). Mathematically, mutually independence of m random variables y_i(t), i = 1, ... m can be expressed by

m    p (y) =  p (y1, ... ,ym) =   Π  p (y i)   i = 1

where p(y) denotes the joint probability density function (pdf) of the random variable y( t ). In other words, signals are independent if their joint pdf can be factorized into marginal distributions.

Second order algorithms are not able to extract or separate random components without temporal structures such as i.i.d. (independent identically distributed) components. Therefore, in this study we used the ThinICA (or Thin SVD ICA) algorithm developed by Cruces and Cichocki (2003). The ThinICA algorithm can be considered as an extension of lower-order algorithms since it employs the second-order and higher-order statistics for the estimation of the rotation matrix U and, consequently, of the demixing matrix W = Â^-1 = U^TQ. The ThinICA algorithm is able to extract simultaneously arbitrary number of components specified by the user. The algorithm is based on a criterion that jointly performs the maximization of several third- and/or fourth-order cumulants of the outputs and/or second-order time-delayed covariance matrices, i.e. on the maximization of the following contrast function:

subject to the constraints U ^T U = I _n , where α _τ are weighting coefficients (typically, equal to 1) and Cum means cumulants for different time tuples t = {t₁, t₂, ... t_p,}. In practice, we have used only the 2-nd, 3-rd and 4-th order cumulants (Cruces and Cichocki, 2003; Cichocki and Amari, 2003).

ICA is a process which statistically reduces a possibly very multidimensional complex data set into sub-components which are statistically independent, and which are expected to capture most of the useful information regarding the underlying events. Since properly estimated ICs are statistically independent from each other, they can be used to create a new set of explanatory variables in order to investigate brain signal relationships more efficiently than it would be possible with the unprocessed data, and can be used by exploratory techniques like DTF or PDC.

3. Results

Following the methodology presented above, we estimated from the high-resolution EEG recordings the cortical waveforms at the ROIs selected for all the group population. The solution of the linear inverse problem returned a series of cortical waveforms for all ROIs considered and for all trials analyzed. In each subject, the cortical waveforms were then subjected to a ThinICA processing, according to the procedures depicted in the Methods section. The independent components obtained by the application of the ThinICA algorithms were then processed by the PDC algorithm, in order to extract the directed causality between the spatial cortical patterns of the estimated data. Such couples of cortical patterns were obtained for each one of the four frequency bands employed in this study.

As an example of the obtained results, Fig. 3 demonstrates the four highest links between the spatial patterns computed for the subject ALFA in the beta frequency band. Each square represents the view of the cortex from the above, nose at the bottom. The causal relationship is expressed between the spatial cortical patterns (squares) located in the same row, so that the left spatial pattern is shown causing the right one. In other words, a spatial cortical pattern related to an independent component located in the left column “causes” the spatial cortical pattern in the right column at the same. The flat-coloured areas in each square depict the cortical segmentation of the ROIs analyzed in this task. The light-purple flat-colored ROI corresponds to parietal cortical area A7, the yellow one – to foot movement area (MIF), the orange ROI - to the primary motor area for the lip movement, the sky blue ROI – to supplementary motor area (SMA), the dark blue ROI – to A6, and the green ROI – to the cingulate motor area (CMA). The colored spheres located in some flat-coloured areas represent the intensity of the cortical activation, as estimated by the ICA (columns of the Anorm matrix). In that way, Fig. 3 presents the strong causal relation in the beta band between the spatial cortical patterns returned by the application of the ICA to the estimated cortical data for subject ALFA. The causality between the spatial cortical patterns in this frequency band is statistically significant at a threshold of 1%.

Note that as the intensity of the spatial pattern of each square is running between 1 and -1, being normalized with the amplitude values moved into the temporal waveforms of the independent components (from the j-th column of the Anorm matrix to the j-th component of the Yscal matrix). Where there are no spheres, no significant cortical activity has been computed. Fig. 4 shows the spatial cortical patterns between the independent components for another normal subject investigated. The issue is now the generation of a series of cortical causality patterns that are “common” between the different subjects.

Figure 3. The figure presents a list of the four strongest links between the spatial patterns computed for the subject ALFA in the beta frequency band. Each square represents the view of the cortex from the above, nose at the bottom. The causal relation is expressed between the spatial cortical patterns (squares) located in the same row with the left spatial pattern causing the right one. In other words, a spatial cortical pattern related to an independent component located in the left column “causes” the spatial cortical pattern in the right column at the same row. This causality is expressed by the blue arrow, moving from the left pattern and pointing to the right one. The flat-coloured areas in each square depict the cortical segmentation of the ROIs analyzed in this task. The sphere located inside each flat-coloured area represents the intensity of the cortical activation, as reflected by the results from the ICA processing (column of the A norm matrix). The figure presents the strong causal relation in the beta band for the subject ALFA between the spatial cortical patterns returned by the application of the ICA to the estimated cortical data. The causality between the spatial cortical patterns in this frequency band is statistically significant at a threshold of 1%. The head below the cortical patterns shows the 3D view corresponding to the square connectivity patterns above.

Figure 4. Same conventions as in Fig. 3 but for the subject ARGI in the beta band.

We applied an appropriate algorithm to locate the similarities between the computed couples of cortical patterns in the different subjects, and the “average” couples of spatial patterns in each frequency band were then generated. Successive figures illustrate these “average” causality patterns due to the preparation for an active foot movement in the different frequency bands analyzed for the group of normal subjects. As already described in the method section, here the yellow spheres indicate that all the subjects in the group have the same ROI engaged in the causality link between cortical patterns, while the red spheres on a ROI indicate that all the subjects but one have the same ROI activated, and the black spheres illustrate that all the subjects but two have the same ROI activated. In the following, only results from cortical causality patterns that presented a spatial correlation of more than 70% are displayed in the different frequency bands.

Fig. 5 shows a cortical pattern in the theta band having 79.9% of correlation between all the subjects.

Figure 5. First average cortical causality pattern in the theta band for a group of normal subjects. Same conventions used as in the previously presented figure. Correlation index between the subjects was 79.9%. The realistic head below shows the equivalent representations on the real cortical surface. The yellow color codes the presence of the ROI in all the subjects of the tested population (6/6), the red color - in 5 out of 6, and the black color - in 4 out 6 subjects.

Taking into consideration the yellow and the red spheres, which indicate the ROI activations in all (yellow), or all but one (red) subject in this study, it seems that in the theta frequency band there was a drive from a cortical network involving the right Cingulate Motor Area (CMA), right foot movement area (MIF), as well as the right parietal area (A7) and supplementary motor area (SMA) toward the right MIF and SMA. Fig. 6 exhibits another independent cortical coupling pattern in the theta frequency band for the analyzed group of subjects.

Figure 6. Second average cortical causality pattern between the independent spatial componenst in the theta frequency band. Correlation between subjects was 72%. Same conventions are used as in the previous figure. Note that the red sphere in the right cyngulate area in the target cortical pattern (on the right) is almost completely hidden from the cortical surface in the realistic reconstruction of the cortex presented here.

In this pattern the MIF bilaterally as well as the right premotor area (A6) were driving the right CMA. No other cortical coupling patterns were statistically significant among the subjects analyzed with more than 70% correlation index in the theta frequency band.

The analysis of the normal subjects exposed a coupling pattern in the alpha frequency band with a correlation index equal to 79%. The pattern is presented in Fig.7.

Figure 7. Causality pattern for the normal subjects in the alpha band. Correlation index was equal to 79.7%.

This pattern revealed that the activation of the primary motor area of the foot bilaterally (MIF), with contributions from the right parietal areas (A7) and CMA, were driving the right MIF and the SMA cortical areas.

We observed two cortical causality patterns over the 70% threshold of the correlation index in the beta frequency bands for the normal subjects in this study. In Fig. 8 the first pattern is shown, which had a correlation index of 71%.

Figure 8. First causality pattern revealed in the normal subject population in the beta frequency band. Correlation index was 71%. Same conventions are used as in the previous figure.

Such a pattern suggests a major involvement of the primary left foot area and the cyngulate motor area in the driving of SMA, bilaterally. The second pattern available for normal subjects in the beta band is characterized by a slightly higher correlation index than before, which means a good agreement between the subjects. In particular, the Figure 9 presents the pattern with a correlation index of 72%.

Figure 9. Second causality pattern revealed in the normalsubject population in the beta frequency band. Correlation index was 72%.

This second causality pattern in the normal group suggests that the primary left foot motor area was driving the activity in the primary motor areas bilaterally.

4. Discussion

The methodological approach presented here attempts to describe the generation of a set of mathematical tools able to depict the existence of distributed cortical networks that drive, or “cause”, in the sense of Granger theory, the activity of other distributed cortical networks, as revealed by analysis of data from high-resolution EEG recordings in humans. This methodological approach uses high-resolution techniques for the estimation of the cortical activity in regions of interest from EEG recordings. The computation of realistic head models for each subject investigated is mandatory in order to generate accurate models of the head and the cortex. By using this technology we were able to estimate the cortical current density in particular ROIs depicted on the realistic model of the cortex and coincident with the Brodmann cortical areas for each subject. Consecutively, the application of ThinICA algorithms for the estimation of independent cortical activations in the ROIs returned a set of time-varying waveforms (temporal independent components) as well as a series of spatial patterns of activations (spatial independent components) that could be further processed. In particular, we analyzed the causal relationship between the independent components in the time domain, by using the established algorithms for the estimation of functional connectivity such as the DTF or PDC to unveil the causality between the ThinICA waveforms. In this way, it was possible to observe patterns of  distributed cortical sources that “caused” in the sense of the Granger theory other patterns of distributed cortical sources in different frequency bands. This means that the inclusion of a particular set of waveforms from the “driving” cortical patterns improved the prediction of the “driven” cortical patterns in the multivariate autoregressive modelling. The presented approach describes how cortical patterns could drive other cortical patterns, by using the concept that such patterns are “independent” in the sense provided by the application of an independent component analysis to the cortical current density estimations. Specifically, the ThinICA method [Cruces-Alvarez and Cichocki, 2003; Cruces-Alvarez et al., 2004; Bakardjian, 2004] was used in order to extract information about the cortical connectivity. In a previous similar approach, Miwakeichi and coworkers (Miwakeichi et al., 2004, 2005) attempted to decompose the space time-frequency components of the EEG during cognitive tasks with the aid of an analysis method called parallel factor analysis (PARAFAC).

Here, the results were provided for a group of normal subjects, performing a real foot movement and the cortical patterns elicited by the above described procedure seems to indicate the existence of a global distributed network in the theta and alpha frequency bands that was driving consistently the right primary motor and supplementary motor areas. In fact, it was possible to observe that the areas activated in all but one subject of the normal population were similar in the theta and alpha band, in both the driving and in the driven networks. The existence of such networks was already underlined by previously existing literature on the simple movements in humans [reviewed in Pfurtscheller and Lopes da Silva, 1999]. However, the simple class of movements investigated here has been chosen because their spatial details were known in advance, thus reducing the uncertainty about the possible results which could be expected. In that way, we were able to test better the proposed methodology that combined independent component analysis (ICA) and functional connectivity estimates.

It is worth of noting that the Granger causality between the cortical areas estimated by the DTF and PDC is not necessarily linked to an increase in the energy or in the amplitude of the potentials of the activated cortical areas. This is due to the fact that the latter observable variables are linked to the synchronized behaviour of the cortical neurons, while the hemodynamic and the metabolic requirements are not [Nunez, 1995; Babiloni et al., 2005]. In fact, the EEG amplitude generated by a patch of cortical tissue, in which only 1% of the neurons discharge synchronously, is 30 times greater than the one produced by the asynchronous discharge of the remaining 99% [Nunez, 1995; Babiloni et al., 2005]. However, the same is not true for the hemodynamic requirements that are rather linked to the firing rate of the same neurons rather than to their firing synchronization.

Another point that should be taken into account, is that the individual spatial components presented in this study are “independent” only in the sense provided by the ThinICA algorithm, while there is no assurance that such networks have a precise physiologic meaning. However, we propose that this assumption could be used as a reasonable working hypothesis in order to enhance the features and the phenomena related to the behaviour and cognition in humans. The methodological approach described in this study could be used in a broader way to investigate not only the connectivity between other cortical areas, but, if faster PDC and DTF algorithms are developed, also to contribute in a novel way to the formation of an advanced Brain-Computer Interface (BCI) with increased number of commands (degrees of freedom) and with improved reliability.

References

Akaike H. (1974) A new look at statistical model identification. IEEE Trans Automat Control AC-19:716-723

Astolfi L., Cincotti F, Babiloni C., Carducci F., Basilisco A, Rossini PM, Salinari S.,. Mattia D., Cerutti S, Ben Dayan D , Ding L., Ni Y, He B. and Babiloni F., Assessing Cortical Functional Connectivity By Linear Inverse Estimation And Directed Transfer Function: Simulations And Application To Real Data, Clin Neurophysiol. 2005 Apr;116(4):920-32

Babiloni F, Babiloni C, Locche L, Cincotti F, Rossini PM, Carducci F. High-resolution electro-encephalogram: source estimates of Laplacian-transformed somatosensory-evoked potentials using a realistic subject head model constructed from magnetic resonance images. Med.Biol.Eng Comput. 2000; 38(5):p. 512-9

Babiloni F., Cincotti F, Babiloni C., Carducci F., Basilisco A, Rossini PM, Mattia D., Astolfi L., Ding L., Ni Y, Cheng K, Christine K, Sweeney J , He B.,  Estimation of the cortical functional connectivity with the multimodal integration of high-resolution EEG and fMRI data by Directed Transfer Function, Neuroimage, 2005 Jan 1;24(1):118-31

Baccalà L.A., Sameshima K, Partial Directed Coherence: a new concept in neural structure determination. Biol Cybern, 84: 463-474, 2001.

Bakardjian H., Common independent components for motion-based brain-computer-interfaces, Proc. of the Society for NeuroScience Meeting, San Diego, USA, Oct. 23-27, 2004.

Brovelli A, Battaglini PP, Naranjo JR, Budai R. Medium-range oscillatory network and the 20-Hz sensorimotor induced potential. Neuroimage. 2002; 16(1):p. 130-41

Buchel C, Friston KJ. Modulation of connectivity in visual pathways by attention: cortical interactions evaluated with structural equation modelling and fMRI. Cereb.Cortex 1997; 7(8):p. 768-78

Cichocki A  and Amari  S, Adaptive Blind Signal and Image Processing: Learning Algorithms and Applications. West Sussex, UK: John Wiley & Sons, 2003.

Clifford Carter G. Coherence and time delay estimation. Proc.I.E.E.E. 1987; 75p. 236-55

Cruces-Alvarez S and Cichocki A, Combining blind source extraction with joint approximate diagonalization: thin algorithms for ICA, in Proceedings of 4th International Symposium on Independent Component Analysis and Blind Signal Separation (ICA2003), (Kyoto, Japan), pp. 463-468, Riken, ICA, Apr. 2003.

Cruces-Alvarez,  S,  Cichocki A, and L. De Lathauwer L, "Thin QR and SVD factorizations for simultaneous blind signal extraction," in Proceedings of the European Signal Processing Conference (EUSIPCO), (Vienna, Austria), pp. 217-220, 2004. (ISBN: 3-200-00165-8).

David O, Cosmelli D, Friston KJ. Evaluation of different measures of functional connectivity using a neural mass model. Neuroimage. 2004; 21(2):p. 659-73

Ding M, Bressler SL, Yang W, Liang H. Short-window spectral analysis of cortical event related potentials by adaptive multivariate autoregressive modeling: data preprocessing, model validation and variability assessment. Biol Cybern 2000;83: 35-45

Franaszczuk P. J., K. J. Blinowska, and M. Kowalczyk, “The application of parametric multichannel spectral estimates in the study of electrical brain activity,” Biol. Cybern., vol. 51, pp. 239–247, 1985.

Friston KJ. Functional and effective connectivity in neuroimaging: a synthesis. Hum.Brain Mapp. 1994; 2p. 56-78

Gerloff C, Richard J, Hadley J, Schulman AE, Honda M, Hallett M. Functional coupling and regional activation of human cortical motor areas during simple, internally paced and externally paced finger movements. Brain 1998; 121 ( Pt 8)p. 1513-31

Gevins AS, Cutillo BA, Bressler SL, Morgan NH, White RM, Illes J, Greer DS. Event-related covariances during a bimanual visuomotor task. II. Preparation and feedback. Electroencephalogr. Clin.Neurophysiol. 1989; 74(2):p. 147-60

Granger C.W.J. Investigating causal relations by econometric models and cross-spectral methods. Econometrica 1969; 37p. 424-8

Grave de Peralta Menendez R, Gonzalez Andino SL. Distributed source models: standard solutions and new developments. In Uhl, C. (ed.). Analysis of neurophysiological brain functioning, Springer Verlag, 1999. pp. 176-201.

He B, Lian J: Spatio-temporal Functional Neuroimaging of Brain Electric Activity, Critical Review of Biomedical Engineering, 30: 283-306, 2002.

Horwitz B. The elusive concept of brain connectivity. Neuroimage. 2003; 19(2 Pt 1):p. 466-70

Inouye T, Iyama A, Shinosaki K, Toi S, Matsumoto Y. Inter-site EEG relationships before widespread epileptiform discharges. Int.J.Neurosci. 1995; 82(1-2):p. 143-53

Jancke L, Loose R, Lutz K, Specht K, Shah NJ. Cortical activations during paced finger-tapping applying visual and auditory pacing stimuli. Brain Res.Cogn Brain Res. 2000; 10(1-2):p. 51-66

Kaminski M, Ding M, Truccolo WA, Bressler SL. Evaluating causal relations in neural systems: granger causality, directed transfer function and statistical assessment of significance. Biol.Cybern. 2001; 85(2):p. 145-57

Kaminski MJ, Blinowska KJ. A new method of the description of the information flow in the brain structures. Biol.Cybern. 1991; 65(3):p. 203-10

Lee L, Harrison LM, Mechelli A. A report of the functional connectivity workshop, Dusseldorf 2002. Neuroimage. 2003; 19(2 Pt 1):p. 457-65

Miwakeichi F, Martinez-Montes E, Valdes-Sosa PA, Nishiyama N, Mizuhara H, Yamaguchi Y, "Decomposing EEG data into space-time-frequency components using parallel factor analysis", Neuroimage 22, pp. 1035-1045 (July 2004).

Mizuhara H, Wang L, Kobayashi K, Yamaguchi Y, "Long-range EEG phase synchronization during an arithmetic task indexes a coherent cortical network simultaneously measured by fMRI", Neuroimage Vol. 27 No. 3, pp. 553-563 (2005)

Nunez, P. L. Neocortical dynamics and human EEG rhythms New York: Oxford University Press. 1995.

Pascual-Marqui RD. Reply to comments by Hamalainen, Ilmoniemi and Nunez. ISBET Newsletter 1995; 6p. 16-28

Pfurtscheller G, Lopes da Silva FH. Event-related EEG/MEG synchronization and desynchronization: basic principles. Clin.Neurophysiol. 1999; 110(11):p. 1842-57

Quian Quiroga R, Kraskov A, Kreuz T, Grassberger P. Performance of different synchronization measures in real data: a case study on electroencephalographic signals. Phys.Rev. 2002; E 65p. 041903

Stam CJ, van Dijk BW. Synchronization likelihood:an unbiased measure of generalized synchronization in multivariate data sets. Physica D 2002; 163p. 236-51

Stam CJ, Breakspear M, van Cappellen van Walsum AM, van Dijk BW. Nonlinear synchronization in EEG and whole-head MEG recordings of healthy subjects. Hum.Brain Mapp. 2003; 19(2):p. 63-78

Taniguchi M, Kato A, Fujita N, Hirata M, Tanaka H, Kihara T, Ninomiya H, Hirabuki N, Nakamura H, Robinson SE, Cheyne D, Yoshimine T. Movement-related desynchronization of the cerebral cortex studied with spatially filtered magnetoencephalography. Neuroimage. 2000; 12(3):p. 298-306

Tononi G, Sporns O, Edelman GM. A measure for brain complexity: relating functional segregation and integration in the nervous system. Proc.Natl.Acad.Sci.U.S.A 1994; 91(11):p. 5033-7

Urbano A, Babiloni C, Onorati P, Babiloni F. Dynamic functional coupling of high-resolution EEG potentials related to unilateral internally triggered one-digit movements. Electroencephalogr.Clin.Neurophysiol. 1998; 106(6):p. 477-87