D ti g r1 g r2 grp

and for multiple time points we can write the compact form:

where p is the number of sources and T the number of time samples. The number n of detectors is typically on the order of 100, while the number of sources p can be either a few (1-10) when using dipole fitting methods or up to several thousand when using cortically constrained imaging methods.

Solutions to the inverse problem can be generally split into two classes: the dipole fitting or scanning solutions and imaging solutions (Darvas, et al., 2004). The aim of inverse methods is typically to minimize the cost function

II I 12

\\D _ G • X 12 with respect to the source parameters, subject to appropriate constraints on the sources X. In the case of dipole fitting, the solution is constrained to consist only of a limited number of sources, which leads to the least-squares dipole fit (Scherg, 1990). Optimizing the subspace correlation of the sources with the data and specifying a correlation threshold instead of the number of sources results in the RAP-MUSIC (recursively applied and projected multiple signal classification) solution of the inverse problem (Mosher, et al., 1999b). Distributed sources over the entire cortical surface or brain volume are allowed in the imaging approaches. To resolve ambiguities resulting from the larger number of sources, the least squares problem is regularized by addition of a norm on the image, which is selected to reflect the power, spatial smoothness or other characteristics of the solution, e.g. Hamalainen, et al., 1993, Pascual-Marqui, et al., 1994. If the regularizer is quadratic then the solution is linear in the data and can be written:

X* = (G • G + C-1) GD, where the matrix C encodes the source constraints. Assuming zero mean noise, the mean of the estimate X can be related to the true source X through the resolution kernel R: E X*J = RX, where R = (0' • G + C-1) 1 G^G. The mean of the estimate of the sources is therefore a linear combination of the true sources.

The beamformer methods (van Veen, et al., 1997, Robinson and Vrba, 1999) represent a hybrid approach to the inverse problem as they make assumptions of temporally independent dipolar sources, but they can also be used to compute activation images throughout the brain volume or on the cortical surface. Alternative solutions to the inverse problem have been proposed such as multipole fits (Jerbi, et al., 2004), which extend the elementary ECD source model to quadrupoles, a minimum norm type solution in a continuous source space (Riera, et al., 1998), and Bayesian and other nonquadratic regularization approaches that emphasize sparseness or other physiologically motivated properties in the solution as reviewed by Baillet, et al. (2001).

2 Statistical Significance of Inverse Solutions from Event Related Data

Once an inverse solution is computed, the question of the statistical significance of the solution arises. In most applications of MEG or EEG, single trial data are contaminated by environmental noise (power lines, electrical equipment), physiological noise (electrocardiogram, eye and muscle artifacts), and spontaneous brain activity. Consequently, meaningful source reconstruction usually requires some form of averaging of the measured event related brain activity over multiple repetitions (or epochs) of the same task. Typically, a pre-stimulus and a post-stimulus time segment is recorded, where the pre-stimulus segment is used to establish a baseline.

The inverse problem is solved either on the event related average of these single epochs or separately for each epoch, depending on the subsequent analysis to be performed. The single epochs, whose number can range from tens to hundreds, can also be used as the basis for assessing the significance of the inverse solution. If the solution is computed by a scanning or dipole fitting method, the significance of individual dipoles can be assessed in terms of the spatial accuracy with which the source is localized (Darvas, et al., 2005, Braun, et al., 1997). If the probability density function of the average data is known, bounds on the spatial confidence intervals for dipoles can be computed analytically by means of the Cramer-Rao lower bounds (Mosher, et al., 1993). Residuals can also be used to assess uncertainty with the chi-squared statistic. However, these methods are highly dependent on the assumed model and do not allow for factors such as missed sources or bias resulting from modeling errors. These problems can be avoided using a nonparametric resampling approach in which we use the single trial data to learn the distribution of the error. One example of such an approach is the Bootstrap method, which can be used to assess dipole localization uncertainty from single trial data (Darvas, et al., 2005). The Bootstrap method approximates the distribution of the data with the sample distribution formed by the collection of single epochs. A new representative sample of the data can then be generated by drawing a new set of single epochs at random and with replacement from the original epochs. Using this technique one can construct confidence intervals for each localized dipole, and reject those for which these intervals do not indicate a sufficiently accurate localization.

In the case of an imaging solution to the inverse problem, significance cannot be readily assessed in terms of localization uncertainty as sources can be reconstructed everywhere in the chosen source space (either cortical surface or brain volume). Instead, we can test at each voxel whether there is significant change in neuronal activity relative to a baseline condition. In more complex experimental designs, significant effects can similarly be assessed by fitting a general linear model at eachvoxel. This form of analysis is very similar to that applied in detecting experimental effects in fMRI data. The problems differ in the degree of spatial correlation in the images and the intrinsic spatial resolution of the two modalities. However, in both cases we need to carefully control for the false positives that may result from performing multiple hypothesis tests (one per voxel). Thresholds for significant voxels in MEG linear minimum norm inverse solutions can be computed using either parameteric random fields methods or nonparametric permutation tests (Pantazis, et al., 2005). In both cases, voxelwise statistics are thresh-olded to control the family wise error rate (FWER), i.e. the probability of one or more false positives. Under the null hypothesis, the maximum distribution of a voxelwise statistic, computed over space or space and time, can be used to select a threshold for a desired FWER. For common distributions (Gaussian, chi-squared, student-t) and sufficiently smooth fields, the upper tail of the maximum distribution can be approximated using the expected value of the Euler characteristic (Worsley, et al., 1996). Application of this approach requires spatial smoothing of the brain activation maps to satisfy random field assumptions. Since resolution in MEG is already rather poor, it is preferable to avoid further smoothing. An alternative nonparameteric approach, which also avoids the need for smoothing, is to use permutations tests to learn the maximum distribution under the null hypothesis. The maximum distribution is learned by randomly permuting the post- and pre-stimulus segments of the epochs. Under the null hypothesis that there is no change in brain activation before and after the stimulus, these should be interchangeable and therefore suitable for generating a new sample of the average dataset. By applying the threshold to the reconstructed image, we create regions of interest that represent statistically significant activity (Fig. 1). Similar non-parametric permutation tests have also been applied to images obtained from beamformers (Singh, et al., 2003).

Time courses of activity for statistically selected regions of interest can be constructed by spatially integrating source activity over the region of interest or selecting the voxel that locally maximizes the test statistic.

In both dipole fitting and imaging solutions, one ends up with a number of localized regions in the brain, each showing significant event related activity. The time series from these locations can then be used for further processing to investigate functional connectivity.

If continuous data are recorded for which no event trigger is available, then the methods as described above cannot be applied. However, steady state activity under different test conditions can be compared using modifications of the thresholding methods described. Alternatively, if anatomical regions of interest are known a priori, then time series can be extracted from inverse solutions at these locations for processing using the methods described below.

Simulated Sources Tikhonov Reconstruction MUSIC Scan LCMV Beamformer Scan

Simulated Sources Tikhonov Reconstruction MUSIC Scan LCMV Beamformer Scan

Fig. 1. Illustration of minimum norm, MUSIC and beamformer inverse methods for two simulated sources and the effect of thresholding on our ability to localize the sources (from Darvas, et al., 2004)

3 Assessing Functional Connectivity from the Inverse Solution

While MEG/EEG does not possess the good spatial resolution of fMRI, these modalities provide an excellent temporal resolution (< 1ms), which fMRI cannot achieve as a result of the relatively slow hemodynamic response function (Kim, et al., 1997). This temporal resolution can be used to look at complex phenomena in the time-frequency domain, such as coherence, phase synchrony, or causality between signals from spatially separated brain regions. The goal of applying such measures is to establish functional correspondence between these regions in order to gain an understanding of the networks of neuronal populations that are involved in executing complex tasks. It should be noted that functional connectivity between regions merely establishes that these regions share mutual information during a specific task, whereas effective connectivity between regions actually implies a causal relationship between these regions (Friston 1994, Lee, et al., 2003, Horwitz, 2003).

Much of the work on functional connectivity using EEG or MEG (e.g. Gevins, et al., 1985, Simoes, et al., 2003) analyzes the connectivity between pairs of electrodes or channels directly. Due to the nature of the electromagnetic forward problem one has to deal carefully with crosstalk of sources. Because the fields or potentials created by an active neuronal source drop off at roughly 1/r2, where r is the distance from the source to the detector, neighboring detectors will typically be sensitive to the same sources, and in some cases all sensors may detect a single source. As described by eq. (2), the signals in each detector will be a weighted linear combination of the signals from all sources. It is possible to reduce sensitivity to crosstalk when analyzing sensor data, either by using connectivity measures that are insensitive to linear crosstalk or through application of statistical tests that allow for crosstalk under the null hypothesis. However, in principle it may be preferable to perform connectivity analysis in the source domain, since the inverse solution can be seen as a spatial unmixing procedure. While the limited resolution of MEG/EEG means that inverse methods will not entirely remove crosstalk between sources, it will certainly be reduced. The maximum amount of crosstalk expected on average from linear inverse methods has been estimated to be less than 40% for EEG data in simulations (Liu, et al., 2002) and will be further reduced by combining MEG and EEG or by including prior information from fMRI. It should be noted that generally the crosstalk will be higher for sources which are located deeper in the brain due to the nonuniform resolution of the inverse solution, as reflected in the resolution kernel R described above.

We now review a number of methods for assessing functional connectivity between brain regions, which are applicable to time series extracted from inverse MEG/EEG solutions.

Covariance

A straightforward approach to investigating dependence between regions is to look at the cross covariance or lagged cross covariance of two signals, the so called Evoked Potential Covariance (EPC) (Gevins and Cutillo, 1993, Smith 1992, Gevins and Bressler, 1988). The EPC is defined as the maximum co-variance over all time lags. While this method can work well on regions that are sufficiently spatially separated, it will be seriously affected by crosstalk and will have very limited spatial resolution. This can be easily demonstrated in the following example:

Let xo(t) and yo(t) be two neuronal sources at locations rx and ry , then the measurements in two channels a(t) and b(t) are given by a(t) = ga (rx) xo (t) + 9a (ry) yo (t), b(t) = gb (rx) xo (t) + gb (ry) yo (t) and the channel covariance is given by cov(a, b) = ga (rx) gb (rx) cov(xo, xo) + ga (ry) gb (r.y) cov(yo, yo)

+ [ga (ry) gb fcj + ga fcj gb (rj] cov(xo, yo)- (3)

Clearly, even if the sources themselves have zero covariance, the channel co-variance will be non-zero, falsely implying a connectivity between the channels. It should be noted, however, that these false interactions can be corrected for in subsequent statistical analysis. Crosstalk can exist regardless of whether there is any true connectivity and therefore will also appear in data under the null hypothesis of no interaction. Consequently a suitable baseline condition can be used to control for false positives. Using a nonparameteric permutation test, it is straightforward to also correct for multiple comparisons (Blair and Karniski 1993).

Coherence

Coherence analysis of two signals is a similar approach in which the signals are represented in the frequency domain. Instead of looking for linear interactions between signals in the time domain, the cross spectrum of the two signals is computed and the coherence found as the magnitude of the cross spectrum normalized by the power spectra of the two signals. Event related coupling in the sensor domain has been widely reported (e.g. Mima, et al.,

2001, Andres and Gerloff, 1999, Miltner, et al., 1999), but similarly to covari-ance analysis, is limited in spatial resolution due to the potential for crosstalk between detectors (Nunez, et al., 1997). Since the cross spectrum can be computed from the covariance of two signals by the Fourier transform, the cross spectrum between two channels can be found by the Fourier transform of eq. (3). Consequently, there is the possibility of falsely inferred interactions, unless controlled for in subsequent statistical analysis. Coherence analysis in the source domain has been used by Gross et al. in 2001 in their method for performing Dynamic Imaging of Coherent Sources (DICS) as the first stage of localizing coherent sources and extracting their time series for further analysis. In DICS the coherence of the sources is computed as the output of a frequency domain beamformer constrained to pass activity with unit gain at a specific location in a specific frequency band. The properties of this frequency domain beamformer are described in detail by Gross, et al. (2003).

The linearly constrained minimum variance (LCMV) beamformer used in DICS has the potential for partial cancellation of coherent sources, particularly if short segments of data are used to find the covariance matrix from which the beamformer weights are computed. Consequently, this method can give inaccurate results for strongly coherent sources. It should be noted, however, that the method is used by Gross et al, (2001) as a first step to localize candidate sources for a subsequent analysis of their phase synchrony (see section phase locking), and that in application of the method to experimental data the actual coherence can be expected to be relatively low (<< 1). Due to crosstalk effects, the original map will usually show a very high coherence value in the vicinity of the reference region and might be misleading. Consequently, when working with cortical coherence maps produced from DICS, it is useful to display the relative increase in coherence as compared to a baseline condition. An example of a typical DICS map with the baseline map subtracted is shown in Fig. 2. Coherence was computed in the 9-14 Hz band for MEG data recorded during a visually cued finger movement of the right index finger.

Since coherence mixes phase and power information, it might not be a suitable measure of connectivity between neuronal assemblies (David, et al.,

2002, Le Van Quyen, et al., 2001, Rodriguez, et al., 1999). Furthermore, the definition of coherence through the cross spectral density requires second order stationarity, which is rarely the case for event related brain signals. We review alternative measures that emphasize phase coupling below.

Fig. 2. Example of a DICS map from 275-channel MEG data for a motor activity experiment in the 9-14 Hz band, where the subject was performing a visually paced movement of the right index finger. The left image shows the coherence map relative to a source in the left Ml region. The right image shows a difference map in which the baseline coherence is subtracted. Both images were thresholded to control the false discovery rate of the DICS output at 1%

Fig. 2. Example of a DICS map from 275-channel MEG data for a motor activity experiment in the 9-14 Hz band, where the subject was performing a visually paced movement of the right index finger. The left image shows the coherence map relative to a source in the left Ml region. The right image shows a difference map in which the baseline coherence is subtracted. Both images were thresholded to control the false discovery rate of the DICS output at 1%

Coherency

Coherency is defined as the imaginary part of the cross spectrum and can be used to detect non-instantaneous interaction between two signals (Nolte et al., 2004). Let x(t) and y(t) be signals computed at two locations in the brain using one of the inverse methods discussed above and assume a set of sources Sj(t) located within the brain. Because of limited resolution, x(t) and y(t) are a linear combination of the true signals. We denote the linear factors as aX and ay. The cross spectrum is then given by:

0 0

Post a comment