Random Fields - Union Intersection tests for detecting time-varying connectivity in EEG/MEG imaging F. Carbonell1 , K. J. Worsley1 , N. J. Trujillo-Barreto2 and R. C. Sotero2

1

Department of Mathematics and Statistics, McGill University, Montreal, Canada 2 Cuban Neuroscience Centre, Havana, Cuba

Abstract Electrophysiological (EEG/MEG) imaging challenges statistics by providing two views of the same spatio-temporal data: Topographic and Tomographic. It is a common practice that Statistical Parametric Mapping (SPM) for these two situations be developed separately. In particular, assessing the statistical significance for spatio-temporal connectivity is a major challenge in this type of studies. This work introduces statistical tests for assessing simultaneously the significance of spatio-temporal correlation structure between ERP/ERF components as well as that of their current generating sources. The multivariate statistical test for detecting connectivity between two sets of scalp EEG/MEG measurements at a given time instant is provided by a greatest root statistic. Consideration of tests at all time instants leads to a multiple comparison problem, which is addressed by the use of Random Field Theory. The Union-Intersection principle is the basis for testing hypotheses about the correlation structure in both topographic and tomographic views. The performance of the proposed method is illustrated with real ERP data obtained from a face recognition experiment.

1

Introduction

Functional connectivity in hemodynamic and electromagnetic data is usually interpreted as the temporal correlation between spatially remote neurophysiological events [9], [10]. This concept can be extended to any other image measurements by saying, for instance, that two different regions of the brain are ”functionally” connected if they show similar anatomical features over subjects [29]. Several approaches have been proposed for assessing functional connectivity in images based on hemodynamic measurements [5], [4], [18], [26], [29], Electrophysiological (EEG/MEG) recordings [12], [25], as well as a multimodal integration data [3]. Typically, these approaches are based on the estimation of the covariance (or correlation) structure between time series measured at different spatial sites. The applicability of any of these approaches is constrained by the obvious limitations of the spatial and temporal resolutions of the measurements. Hence, those recording modalities producing image data with high spatial and temporal resolution are very suitable for an effective dynamic (time-varying) functional connectivity analysis. Recent years have seen the emergence of a number of methods for high resolution EEG/MEG data, characterized by the use of a large amount of scalp electrodes. They hold the promise of complementing other imaging modalities, such as PET or fMRI, by providing very high temporal resolution maps of neural activation. These techniques however challenge the statistics of functional connectivity in several ways. One problem is the need to detect significant correlations between the 1

components of the Event Related Potential/Event Related Field (ERP/ERF) and their topographic distribution, that is, the spatio-temporal correlation structure of the ERP/ERF components in high dimensional EEG/MEG data. In addition, interest is increasingly focusing on the identification of the current sources that generate such components. As a consequence, Tomographic views of the same data resulting from the solution of a linear inverse problem [21] require also of a much more complex connectivity analysis. Different approaches have been followed for the analysis of connectivity in both, topographic and the tomographic views of EEG/MEG data. A number of linear and nonlinear connectivity measures have been proposed for analyzing the spatio-temporal topographic correlation structure of ERP/ERF components (see e.g. [8] and [19]). Examples are multichannel EEG/MEG crosscorrelation coefficients (or coherence in the frequency domain) [11], [22], Directed Transfer Function (DFT) [16], mutual information, nonlinear correlation, etc. On the other hand, some of these measures has been adapted for analyzing connectivity over the current density sources obtained by solving an EEG/MEG inverse problem [1], [23]. Indeed, for the connectivity analysis in the tomographic view one could use, in principle, any of the available methods for PET, MRI and fMRI images. The main problem with these approaches is that they use different criteria for topography and for tomography, in fact employing statistical tests derived from different principles. Besides, almost none of these methods provide a decision threshold for statistical significance in a classical hypothesis testing framework. As a consequence there is no general method for assessing simultaneously the spatio-temporal correlation structure of both, the ERP/ERF components and theirs generating current density sources. This work presents a blend of techniques that provides a global decision threshold for detecting changes in the time-dependent correlation structure between ERP/ERF components recorded on two different scalp regions. It is shown that the same threshold can be used for assessing the statistical significance of both the topographic and tomographic correlation structure of such components. In particular we are going to apply the classical Statistical Parametric Mapping (SPM) approach based on Random Field Theory (RFT). As it was pointed out in [17], there exists a severe limitation when using space-time SPM for ERP/ERF studies. For instance, it is very difficult to make inferences about the temporal extent of the evoked responses within the spatiotemporal SPM framework. So, in most of the cases it is necessary inferring over time explicitly. That is just the approach followed here. We model our spatio-temporal ERP/ERF data as a multivariate response at each time instant and make inferences in the context of ”temporal” SPM. Specifically, the statistical procedure proposed in this paper can be summarized in the following steps: 1)Definition of the global null hypothesis of ”no temporal connectivity” between the ERP/ERF components generated in two different brain regions. 2)Calculation of a multivariate greatest root statistic θt at each time instant t in the time window of analysis for testing the hypothesis of Step 1. 3)Computation, by RFT, of a global decision threshold Gα (with the assumption that the statistics of Step 2 as realizations of a 1-dimensional random field). Detection of ”connectivity” between ERP/ERF components at those time instants for which θt exceeds Gα . 4)Application of the Union-Intersection (UI) tests for examining hypotheses about the topographic and tomographic correlation structure (at the time instants selected by Step 3) of the ERP/ERF components. This main idea underlying this approach is that, while EEG/MEG data is usually represented by vectors of voltage/field values recorded from a set of sensors, many Topographic and Tomographic 2

views of the same data are obtained by linear combinations (LC) of these vectors. For example, different montages of the EEG such as bipolar or Laplacian are just LCs of the original data. Additionally, Linear Inverse Solutions are also LCs of the original data. The UI principle provides tests for any chosen LC. Thus, the power of this proposal is that one could construct statistical tests for the topographic correlation structure of the original sensor recordings. On the other hand, it is also possible to provide statistical tests for the correlation structure in all LC of generators produced by any conceivable linear EEG/MEG inverse solution. The use of UI principle is not new in the context of statistical analysis of neuroimaging data. In fact, the statistical methodology presented in this paper is an extension of the analysis carried out in [6]. In that paper, RFT and UI principle were used almost in the same way as the approach followed here. There, RFT solved the multiple comparisons problem resulting from testing the global null hypothesis of ”no ERP/ERF components” , while UI tests provided the tools for assessing the topographic and tomographic distribution of the detected ERP/ERF components. In this way, the method followed in the present paper takes an insight to the analysis of the spatial-temporal correlation structure of ERP/ERF components in EEG/MEG measurements. More recently, UI principle has been also used in the problem of finding the P-value of the maximum of random field of some multivariate test statistics, with applications in brain shape analysis [24].

2

Methods

Let [X]sct be a three-dimensional array that represents the recorded voltages values corresponding to trial s, s = 1...S, channel c, c = 1...C and time instant t, t = 1...T. Here, S, C and T denote, the number of trials, channels and time instants, respectively. The common model for such ERP data [2] assumes that the recorded voltage values [X]sct are the superposition of the ERP (“signal” µct ) and a background brain activity “noise” (esct ) noncontingent to the stimuli, that is, Xsct = µct + esct , where the vectors est = (es1t , ..., esCt ) are assumed statistically independent across trials and normally distributed with mean equal to zero and (typically unknown) covariance matrix Σt .

2.1

Statistics for global time-dependent connectivity

As in the classical problem of detecting temporal ERP/ERF components, it would be interesting to give an analogous answer for the problem of detecting ”temporal” connectivity between such components. That is, detecting changes in the time-dependent spatial correlation structure of ERP/ERF components. For this problem, one would like to test, in principle, the null hypothesis that the covariance matrix Σt is diagonal for all time instants t. However, it is well-known that the relationship between current source densities localized in the brain and the scalp electric field recorded in EEG/MEG measurements is determined by the volume conduction properties of the tissues in the head and its properties. Hence, such global null hypothesis for Σt being a diagonal matrix would be almost always rejected, mainly in the case of high density multichannel recordings. That is, even the single case of a bipolar source inside the brain would produce an electrical field extended over a large region on the scalp surface and as a consequence, resulting in high correlation values between spatially-closed electrodes located in such region.

3

Therefore, a more realistic hypothesis testing to deal with would be a one that involves timedependent connectivity between ”spatially remote” regions, for instance, inter-hemispheric connectivity at each time instant. The mathematical formulation for this is as follows. Partition the C channels into two disjoint sets I1 and I2 with cardinality C1 and C2 , respectively, where C1 +C2 = C and I1 ∪ I2 = set of all channels. In this way, the covariance matrix Σt is partitioned as well as µ 11 ¶ Σt Σ12 t Σt = , (1) Σ21 Σ22 t t 22 where Σ11 t and Σt are square matrices of dimension C1 and C2 , respectively. In order to detect connectivity between the channels of each subset of the partition, one would like to test the global null hypothesis

H0 : Σ12 t = 0, for all t = 1...T, which in turn can be decomposed in the following T marginal null hypotheses: H0t : Σ12 t = 0, for each t = 1...T. The likelihood ratio test for H0t is a Wilks’s Lambda statistic [20] given by Λt = where

µ St =

21 11 −1 12 det[(S22 t − St (St ) St )] , det[(S22 t )]

S11 S12 t t 21 St S22 t

¶

S 1X (xst − xt )(xst − xt )0 = S s=1

is the maximum likelihood estimate of the covariance matrix Σt , with xst = (Xs1t , ...XsCt )0 and xt =

S 1X xst . S s=1

According to ([20], pp. 135), under H0t , Λt follows a Wilk’s Lambda probability distribution Λ(C2 , S − 1 − C1 , C1 ).

2.2

Random Field Theory

For a fixed time instant t, a test for H0t at the significance level α is given by rejecting this hypothesis when Λt is greater than the (1 − α)100 percentile of the distribution Λ(C2 , S − 1 − C1 , C1 ), say cα . However, a test for H0 would require the application of several tests, each one at the significance level α, which leads to a well known multiple comparisons problem [15]. The classical strategy to solve this problem is to use a Bonferroni-type correction. However, this is usually an excessively conservative procedure and leads to an unnecessary loss of statistical power for the case of correlated tests, which is likely to occur in our setting due to the temporal correlation structure that usually appears in multi-channel EEG/MEG data. A more effective solution to the correlated multiple comparisons problem is the use of RFT for generating distributional approximations for the maximum of random field statistics. The application of this approach to the current setting reduces to consider that the Λt statistics are realizations from a (1D) Wilk’s Lambda random field over the interval(search region) U . Unfortunately, the random field theory for this statistics has not so far fully studied. 4

However, the Roy’s maximum root statistic has been reported to be a successful alternative to the Wilk’s Lambda statistics [28], [24]. The Roy’s maximum root random field λt , t ∈ U , defined as in h 22.1 i−1 h 22 22.1 i Mt Mt −Mt [24], is given by the largest eigenvalue of the matrix S−1−C , where C1 1 21 11 −1 12 M22.1 = S(S22 t t − St (St ) St ), M22 = S(S22 t t ).

However, we used instead the so-called greatest root statistic θt , (defined on page 84 in [20]), which is related to the Roy’s maximum root statistic by θt =

et λ et ) (1 + λ

,

et = C1 λt . Our choice is based on the fact that the greatest root statistic θt is just but where λ S−1−C1 the Union-Intersection test statistic associated to the null hypothesis H0t (page 136 in [20]). As we shall see in the next subsection, this fact has useful consequences for the rest of the paper. For the −1 21 11 −1 12 meantime, note that θt is the largest eigenvalue of the matrix (S22 t ) St (St ) St , this is, −1 21 11 −1 12 θt := max eig((S22 t ) St (St ) St ),

which follows a probability distribution denoted by θ(C2 , S − 1 − C1 , C1 ). Therefore, a test for H0 shall be based on the statistic θmax := max(θ1 , ..., θT ), where H0 is rejected when θmax > Gα . Here, Gα denotes the (1 − α)100 percentile of the distribution of θmax under H0 , which can be calculated according to the results presented in [24], which, for completeness we have reproduced in the Appendix. Finally, the connectivity revealed during the generation of the ERP/ERF components are detected by the conex components of the excursion set of the random field θmax above the threshold Gα . This is, at those intervals [t1 , t2 ] ⊂ U for which θt > Gα holds for all t ∈ [t1 , t2 ].

2.3

Connectivity on topographic view

According to the previous subsection, those time instants t for which θt > Gα are the more interesting ones. Let τ be one of those time instants (i.e. θτ > Gα ). One of the practical consequences of the Union-Intersection (UI) tests is that, if H0 is rejected then one can ask which pairs of entries in the matrix Σ12 τ were responsible ([20], pp. 136). For this purpose the following hypotheses are considered: ab H0τi j : (Σ12 τ )ij = 0, where the vectors ai = (0, 0, ....1, ...0), bj = (0, 0, ....1, ...0), i = 1, ..., C1 , j = 1, ..., C2 have a number 1 in the i − th and j − th components, respectively. The use of the above notation comes from the following fact. According to the given partition, the partitioned array data X1 = [X]s,c1 ,τ , c1 ∈ I1 and X2 = [X]s,c2 ,τ , c2 ∈ I2 , s = 1, ..., S, can be considered as sample data matrices of multidimensional vectors x1 and x2 , respectively. Hence, testing (Σ12 τ )ij = 0 is equivalent to test that the sample correlation coefficient rai bj between the vectors a0i x1 and b0j x2 is zero. This is, ra2i bj =

2 [C(a0i x1 , b0j x2 )]2 [a0i S12 τ bj ] = 0 22 V (a0i x1 )V (b0j x2 ) [a0i S11 τ ai ][bj Sτ bj ]

5

2 Hence, the UI test based on this decomposition uses the statistics max{ra,b }, which equals the a,b

−1 21 11 −1 12 largest eigenvalue of (S22 τ ) Sτ (Sτ ) Sτ ([20], pp. 136), this is, the greatest maximum root ab statistics with distribution θ(C2 , S − 1 − C1 , C1 ). Therefore, according to the UI principle, H0τi j is rejected when ra2i bj ≥ Gα .

2.4

Connectivity on tomographic view

As it is well known, one should be cautious when interpreting the ERP topography in terms of neural activation regions, which requires the use of specific tomographic techniques. Indeed, the development of physical models and mathematical methods has led to a novel imaging technique, usually called Brain Electromagnetic Tomography (BET), which has become a very useful tool for determining the possible neural generators of ERP components. It is based on the estimation of the Primary Current Density (PCD) Jτ , given the EEG xτ recorded on the scalp surface (in our seeting xτ represents a multidimensional vector of the sample data matrix [X]s,c,τ , s = 1, ..., S, c = 1, ..., C). This is known as the inverse problem of the EEG and requires the solution of the system of linear equations xτ = KJτ + eτ , where K is known as the leadfield and eτ is a vector of measurement errors. Although the methods developed in this paper are valid for any linear inverse solution, their use shall be illustrated with the popular LORETA (Low Resolution Electromagnetic Tomography) solution [21]. Under Gaussianity and independency assumptions for the noise eτ , and smoothness constraints for Jτ (in the sense of the second derivative), the LORETA solution can be expressed as b τ ) = (K0 K + r2 L0 L)−1 Kxτ , J(x b τ ) is the estimated PCD, r is a regularizing parameter and L is a discrete version of the where J(x Laplacian operator. In order to clarify the exposition, in the sequel we shall omit the dependence on τ . Suppose b b1 (x)0 ...., J bng (x)0 ) has dimension 3ng , where ng denotes the number of that the vector J(x) = (J bg (x), g = 1 . . . ng represents the three generators sampled over a 3-D domain in the brain and J spatial components of the current density for the g − th generator. From the expression given for b bg (x) is a certain linear transformation Ag x of the vector J(x), it is easily seen that each vector J bg (x). x. Therefore, hypotheses about Ag x are equivalent to the corresponding hypotheses about J In principle, we are interested on testing for null correlation between any pair of generators g and bg (x) and J bh (x)). However, our analysis shall be restricted to h (i.e. all correlations between J corresponding generators obtained by the original partition given in vector x. That is, we will be doing a tomographic connectivity analysis for seeing the influence of the generating sources in each subset of channels I1 and I2 . Mathematically, this means that we are going to be testing for null bg (x1 ) and J bh (x2 ) , g, h = 1, ..., ng . For g and h fixed, this correlations between all possible pairs J task involves the sample correlation between all the possible pair of components of the 3-D vectors A1g x1 and A2h x2 (for certain matrices A1g and A2h ), where, as in the previous subsection, x1 and x2 are multidimensional vectors. This is the typical statistical problem of canonical correlation. So, we are be testing that the sample maximal canonical correlation between the vectors A1g x1 and A2h x2 be zero, which is defined as 2 rA 1 ,A2 = max g h

u,v

[C(u0 A1g x1 , v0 A2h x2 )]2 . V (u0 A1g x1 )V (v0 A2h x2 ) 6

Notice that rA1g ,A2h equals the maximal eigenvalue of the matrix (Y20 Y2 )−1 Y20 Y1 (Y10 Y1 )−1 Y10 Y2 , where Y1 0 = A1g X1 0 and Y2 0 = A1h X1 0 , with Xl = [X]s,cl ,τ , cl ∈ Il , l = 1, 2. Hence, as in the 2 previous subsection, UI tests are based on the statistics max {rA 1 ,A2 }, which equals the largest 1 2 A ,A

−1 21 11 −1 12 eigenvalue of (S22 τ ) Sτ (Sτ ) Sτ . Thus, according to the UI principle, the null hypothesis of zero bg (x1 ) and J bh (x2 ) is rejected when r2 1 2 ≥ Gα . maximal canonical correlation between J Ag ,Ah From a practical point of view, significant connectivity is obtained by thresholding at Gα the bg (x1 ) and J bh (x2 ) , g, h = 1, ..., ng . However, 6D random field of all cross-correlations between J it requires a lot of computational effort and a carefully interpretation of the results. This comes from the fact that, likely, the significant cross-correlations be those resulting from common sources generated in both subset of channels. Therefore, our results in the next section are restricted to the analysis of just the 3D ”homologous ” correlation random field [5] obtained by the correlations bg (x1 ) and J bg (x2 ) , g = 1, ..., ng . between J

3

Applications

As an illustration of the developed method, we apply our results to an ERP data set of a face recognition paradigm used in [14]. As it was described in [14], the basic paradigm consists on presentation of 86 faces (43 famous and 43 novels) and 86 scrambled faces, which gives 3 event-types: unfamiliar faces (U), familiar faces (F) and scrambled faces (S). The scrambled faces were created by 2D fourier transformation, random phase permutation, inverse transformation and outline-masking of each face. Thus faces and scrambled faces are closely matched for low-level visual properties such spatial frequency power density. Faces were presented for 600ms, every 3600ms. More details about the experimental paradigm can be found at www.fil.ion.ucl.ac.uk/spm/data/mmfaces.html. The EEG data were acquired on a 128-channel Biosemi ActiveTwo system, sampled at 2048 Hz, plus electrodes on left earlobe, right earlobe, and two each to measure HEOG and VEOG. The data were referenced to the average of left and right earlobe electrodes and epoched from -200ms to +600ms. These epochs were then detrended and examined for artifacts, defined as time points that exceeded an absolute threshold of 120 microvolt (mainly in the VEOG). After artifacts rejection, our primary data set on conditions F, U and S consisted on 38, 38 and 76 trials of 128 channels each. To get an idea of the data, the overall mean (over all 76 differences face-scrambled face) of the EEG data at 8 of the 128 electrodes is shown in Figure 1(a). As it was pointed out in [14], it can be seen some clear differences between faces (F and U) and scrambled faces (S) that appear as a negative ERP component around 110 ms for the electrodes CZ, PZ, C4, and FPZ, FZ, located in the temporal and frontal region, respectively. Note also that this negative component is followed by an enhancement of a positive peak around 170 ms. It is also evident an enhancement of a negative component around 170 ms at occipito-temporal channels IZ, P8 and P7 (the so-called N170 component), which in turn is preceded by a positive peak around 110 ms. These facts can be also seen in the topographic maps showed in Figure 1(b). As it was described in Section 2, our statistical tests assume the decomposition of the scalp channels set into two disjoint subsets. That is, the raw data is partitioned in two array data [X]s,c1 ,τ , c1 ∈ I1 and [X]s,c2 ,τ , c2 ∈ I2 , where in this example, I1 and I2 are sets of scalp channels corresponding to the left and right hemispheres, respectively. To remove the effect of familiar and unfamiliar faces depicted in Figure 1, we used the residuals obtained by subtracting the average of each group (familiar, unfamiliar) from the data in that group, to give S = 76 − 2 = 74 effectively independent observations. 7

Figure 1: (a) Left: Average (over 76 trials) of the difference face - scrambled face EEG data for the channels CZ, PZ, IZ, P8, C4, FPZ, FZ and P7, plotted against time t(ms). There is a negative EEG peak around t = 110ms for the electrodes CZ, PZ, C4 in the temporal region, and FPZ and FZ in the frontal region, followed by a positive peak around 170ms. The reverse happens in the occipito-temporal channels OZ, P8 and P7 (the so-called N170 component). (b) Right: the scalp topographic distribution of the 128 channels at the time instants t = 110ms and t = 170ms. Figure 2(a) shows the resulting greatest root statistic θt for t ∈ U = [−200, 600] ms and the decision threshold (dash vertical line) obtained by random field theory as described in the previous section. For this analysis, just 11 scalp channels (those shown in Figure 2(b)) per hemispheres were selected (see a brief discussion about this choice in Section 4). Notice the presence of several peaks of significative local maximum greatest root statistics. Our subsequent analysis focus on the time intervals containing some significative instants marked with circles in the plot. As expected, there is significative global connectivity between both hemispheres in the time interval [100, 190]ms. However, some other interesting time intervals with significative values of greatest root statistics are [280, 330]ms, [390, 450]ms and [480, 510]ms. These finding seem to be related with the activity produced by the N170 component and certain significant activation onsetting 300ms, usually associated to a familiarity effect in the face recognition process [14]. The analysis of the topographic connectivity maps reveals no significant correlation at the time instant t = 111ms and a significant correlation at t = 170ms, which should be associated to the bilateral occipital maximal negativity observed in Figure 1 (b). This result agrees with those of [14], where it was reported a bilateral activation in the occipital region around t = 170ms, associated 8

Figure 2: (a)Top: Greatest root statistics and decision threshold for an inter-hemispheric comparison. (b)Bottom: Topographic maps of the correlation analysis at four different time instants of local maximum greatest maximum root statistics (marked with circle in the top plot) to a face perception process. Notice also a bilateral connectivity between occipital and temporal regions at t = 322ms and t = 427ms, as previously commented, related to certain familiarity effect. Tomographic connectivity maps at these four time instants are displayed in Figure 4. Figures 4a, 4b, 4c and 4d corresponds to the time instants 111ms, 170ms, 322ms and 427ms, respectively. Notice that there exists a consistent significant homologous correlation at the right fronto-temporal region, although this is extended to the right occipital region at time instant 170 ms. This fact suggests that such regions are involved in the face recognition process as generated sources of the selected scalp measurements in both hemispheres. This conclusions agrees with that reported in [14], where it was suggested that the brain activation produced by the face recognition process is mainly right lateralized. These results are also consistent with the classical cognitive ”core model” for face recognition and perception [13], since this cognitive process involves the processing and transfer of information in a neuronal network that includes brain structures in the occipital and frontal regions as well as in the superior-temporal area.

9

(a)

(b)

(c)

(d)

Figure 3: Tomographic maps of the homologous correlation analysis between the sources generated by the left and right scalp channels at times t = (a) 111ms, (b) 170ms, (c) 322ms, (d) 427ms. Letters A, P, L, R in each subplot mean Anterior, Posterior, Left and Right, respectively. We conclude that there are statistically significant homologous correlations mainly in the right fronto-temporal region

10

4

Discussion

In this section we briefly present some extensions and some limitations of our proposed method. As one realize, the exposition of our method is based on a connectivity analysis within a single EEG/MEG trial(condition) type. However, the whole analysis can be easily extended for comparison between two trials type. Indeed, in this former case, the array data [X]s,c1 ,t , c1 ∈ I1 and [X]s,c2 ,t , c2 ∈ I2 would be considered as multivariate sample data matrices on two different trials type. Even more, they could be also considered as ERP/ERF data obtained from subjects on two different experimental conditions. These kind of comparisons would lead to more meaningful applications and interpretation of the results when dealing with different experimental conditions (contrasts). A possible limitation of the methods presented here is the basic model for the ERP/ERF as the linear superposition of the event related activity and background activity. This model is only a first approximation, useful for exploratory purposes. More realistic, nonlinear and non-stationary models of the ERP/ERF may require more complex statistical procedures. The extension of the test statistic to a whole time interval hinges on the use of RFT. As it was pointed out in [24], 2ρd (λ) is not the EC density of a Roy’s maximum root random field, rather it is twice the alternating sum of the EC density of all roots. But, the authors in ([24]) consider that, for high thresholds, the other roots are much more smaller than the maximum and so, theirs EC is close to zero. This is the reason for the 1/2 that appears in formula (2). Nevertheless, some simulations carried out in the preparation of this paper suggest that this last assumption is quite strong for Roy’s maximum random fields of high dimensionality. In this paper, the dimension in the greatest maximum root distribution θ(C2 , S − 1 − C1 , C1 ) is given by C2 , where without loss of generality we have considered C2 ≤ C1 . So, for high values of C2 (values greater than 15) our simulations turned out in high values (too conservative) of the decision threshold Gα as well (to our knowledge, the Roy’s maximum root random field results have just been applied in brain shape analysis with C2 = 3). These considerations leaded us to choose the value C2 = 11 in Section 3 (i.e. just 11 scalp channels per hemisphere). In principle one could consider higher dimensionality by doing the whole (128 channels)inter-hemispheric connectivity analysis with C2 ≈ 60. However, we anticipate that such type of analysis would also leads to extremely high values of the greatest root statistics over the whole time window. This indicate that our greatest root statistics is perhaps sensitive to the extant high correlations that appear between scalp channels spatially closed. Hence, a more meaningful analysis in the topographic and tomographic view would be required. In such situation we then recommend the use of a recent result for time-varying cross-correlation (2D or 3D) random fields [7].

5

Conclusions

A novel application of the RFT and the UI tests was introduced in the framework of EEG/MEG hypothesis-testing problems. RFT was applied for solving the multiple comparisons problem that appears when examining the global null hypothesis of “no temporal connectivity”. On the other hand, UI tests were used for testing hypothesis about the topographic and tomographic distribution of the detected temporal connectivity between ERP/ERF components. Therefore, the combination of both techniques provided a general criterion for assessing simultaneously the spatio-temporal detection of the correlation structure between ERP/ERF components as well as the time dependent cross-correlation between its generating current density sources. Finally, physiologically meaningful results were obtained in the application of the statistical methodology to actual EEG data recorded from a face recognition experiment. 11

6

Appendix

Computation of Gα For any value θ, the explicit approximation for P (θmax ≥ θ) is given by P (θmax ≥ θ) ≈ ρ0 (λ) + Resels1 ρ1 (λ), ³ where λ =

S−1−1C1 C1

´¡

θ 1−θ

¢

and Roy’s maximum EC densities ρd (λ), d = 0, 1 are given by (see [24]) C

2 1X ρd (λ) = µj (UC2 )ρFd+j (λ). 2 j=0 C +1

2 2j+1 π j/2 Γ( 2 ) C +1−j j! Γ( 2 2 )

In the formula above, µj (UC2 ) = µ ρFd (λ) =

log 2 π

¶ d2

min(i,d−1−i)

×

X

and the F EC densities are given by

2(d − 1)!Γ( p+m−d ) 2 m

p−d 2

µ

Γ( p2 )Γ( m2 )

pλ 1+ m

¡ p+m−d +j−1¢¡m−1¢¡ 2

j

(2)

i−j

¶− p+m−2 d−1 X 2 p−d (−1)d−1−i (pλ)i+ 2

p−1 d−1−i−j

i=0

¢

m−i ,

j=0

with p = C1 and m = S − 1 − C1 degrees of freedom. The number of 1D resels, Resels1 , can be estimated unbaisedly [27] as follows. Let rct be the sample correlation coefficient between Xs,c,t and Xs,c,t−1 over s for fixed c, t. Then Resels1 = (4 log 2)

− 21

C T 1 XXp 2(1 − rct ). C c=1 t=2

Finally, the threshold tα is obtained by solving the equation P (θmax ≥ Gα ) = α.

Acknowledgement The authors would like to thank Rik Henson and Will Penny for kindly providing the experimental data.

References [1] Astolfi, L., Cincotti, F. Mattia, D., Marciani, G. Baccala, L. Fallani, F., Salinari, S., Ursino, M., Zavaglia, M., Ding, L. Edgar, J. C., Miller, G. A., He, B. & Babiloni, F. (2007). Comparison of Different Cortical Connectivity Estimators for High-Resolution EEG Recordings. Human Brain Mapping 28,143-157. [2] Aunon, J.I., McGillen, C.D. & Childers, D.C. (1981). Signal processing in evoked potentials research: averaging and modeling. CRC Critical Reviews in Bioengineering 5, 323-367.

12

[3] Babiloni, F., Cincotti, F., Babiloni, C., Carducci, F., Mattia, D., Astolfi, L., Basilisco, A., Rossini, P. M., Ding, L., Ni, Y., Cheng, J., Christine, K., Sweeney, J. & Heg, B. (2005). Estimation of the cortical functional connectivity with the multimodal integration of highresolution EEG and fMRI data by directed transfer function. NeuroImage 24, 118-131. [4] Baumgartner, R., Ryner, L., Richter, W., Summers, R., Jarmasz, M. & Somorjai, R. (2000). Comparison of two exploratory data analysis methods for fMRI: fuzzy clustering vs. principal component analysis. Magnetic Resonance Imaging 18, 89-94. [5] Cao, J. & Worsley, K.J. (1999). The geometry of correlation fields with an application to functional connectivity of the brain. Annals of Applied Probability 9, 1021-1057. [6] Carbonell, F., Gal´an, L., Vald´es, P., Worsley, K. J., Biscay, R. J., D´ıaz-Comas, L., Bobes, M. A. & Parra, M. (2004). Random Field-Union Intersection Tests for Linear Electrophysiological Imaging, Neuroimage 22, 268-276. [7] Carbonell, F., Worsley, K. J., Trujillo-Barreto, N. & Vega, M. (2007). The geometry of timevarying cross-correlation random field, submmited. [8] David, O., Cosmelli, D. & Friston, K. J. (2004). Evaluation of different measures of functional connectivity using a neural mass model. NeuroImage 21, 659-673. [9] Friston, K.J., Frith, C.D., Liddle, P.F., & Frackowiak, R.S. (1993). Functional connectivity: the principal-component analysis of large (PET) data sets. J. Cereb. Blood Flow Metab. 13 (1), 5-14. [10] Friston, K. J. (1996). Statistical parametric mapping and other analyses of functional imaging data. In: Toga AW, Mazziotta JC, editors. Brain mapping: the methods. San Diego: Academic Press, 363-396. [11] Gevins, A. & Remond, A. (1987). Methods of analysis of brain electrical and magnetic signals. Handbook of Electroencephalography and Clinical Neurophysiology, Vol1, Amsterdam, Elsevier. [12] Gevins, A. (1989). Dynamic functional topography of cognitive task. Brain Topography 2, 37-56. [13] Gobbini, M.I. & Haxby, J.V. (2007). Neural systems for recognition of familiar faces. Neuropsychologia 45(1), 32-44. [14] Henson, R.N., Goshen-Gottstein, Y., Ganel, T., Otten, L. J., Quayle, A. & Rugg, M. D. (2003). Electrophysiological and Haemodynamic Correlates of Face Perception, Recognition and Priming. Cerebral Cortex 13, 793-805. [15] Hochberg, Y. & Tamhane, A. C. (1987). Multiple comparisons procedures. New York: Wiley. [16] Kaminski, M., Ding, M., Truccolo, W. A. & Bressler, S. (2001). Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of significance. Biological Cybernetics 85, 145-157. [17] Stefan J. Kiebel & Karl J. Friston. (2004). Statistical parametric mapping for event-related potentials: I. Generic considerations. NeuroImage 22, 492-502. 13

[18] Koch, M.A., Norris, D.G., & Hund-Georgiadis, M. (2002). An investigation of functional and anatomical connectivity using magnetic resonance imaging. NeuroImage 16, 241-250. [19] Koenig, T., Studer, D., Hubl, D., Melie, L. & Strik, W. K. (2005). Brain connectivity at different time-scales measured with EEG. Phil. Trans. R. Soc. B 360, 1015-1024. [20] Mardia, K. V., Kent, J. T. & Bibby, J. M. (1979). Multivariate Analysis. Academic Press. [21] Pascual-Marqui, R. D., Esslen, M., Kochi, K. & Lehmann, D. (2002). Functional imaging with low-resolution brain electromagnetic tomography (LORETA): a review. Methods and Findings in Experimental and Clinical Pharmacology 24C, 91-95. [22] Kaspar Schindler, Howan Leung, Christian E. Elger & Klaus Lehnertz. (2007). Assessing seizure dynamics by analysing the correlation structure of multichannel intracranial EEG. Brain 130, 65-77. [23] Thatcher, R. W., Biver, C. J. & North, D. (2007). Spatial-temporal current source correlation and cortical connectivity. Clinical EEG Neuroscience 38(1), 35-48. [24] Taylor, J.E. & Worsley, K.J. (2007). Random fields of multivariate test statistics, with applications to shape analysis and fMRI. Annals of Statistics, accepted. [25] Urbano, A., Babiloni, C., Onorati, P. & Babiloni, F. (1998). Dynamic functional coupling of high resolution EEG potentials related to unilateral internally triggered one-digit movements. Electroencephalogr. Clin. Neurophysiol. 106 (6), 477-487. [26] van de Ven, V.G., Formisano, E., Prvulovic, D., Roeder, C.H. & Linden, D.E. (2004). Functional connectivity as revealed by spatial independent component analysis of fMRI measurements during rest. Human Brain Mapping 22, 165-178. [27] Worsley, K.J., Andermann, M., Koulis, T., MacDonald, D. & Evans, A.C. (1999). Detecting changes in non-isotropic images. Human Brain Mapping 8, 98-101. [28] Worsley, K.J., Taylor, J.E., Tomaiuolo, F. & Lerch, J. (2004). Unified univariate and multivariate random field theory. NeuroImage 23, S189-195. [29] Worsley, K.J., Chen, J-I., Lerch, J. & Evans, A.C. (2005). Comparing connectivity via thresholding correlations and SVD. Philosophical Transactions of the Royal Society 360, 913-920.

14