T

Phase Locking

An alternative to the classical coherence measure is to separate phase and amplitude information in the signals and analyze the phase locking or phase synchrony as a measure of connectivity between neural assemblies. Phase locking between two oscillators is defined as \^n,m(t)\ < const,pn,m(t) = n^i(t) — mfi2(t), where n and m are integers and ^i(t) and (t) are the instantaneous phases of the signals (Tass, et al., 1998). In order to compute any measure of phase synchrony, the instantaneous phase of the signals must first be computed. This can be done by first narrowband filtering of the signals and then extracting the instantaneous phase by applying the Hilbert transform (Tass, et al., 1998). Alternatively phase can be computed from a time frequency decomposition using Gabor wavelets (Lachaux, et al., 1999). Once the instantaneous phase has been extracted from the signals, the n:m synchronization index (Tass, et al., 1998) or the phase locking value (PLV) (Lachaux, et al., 1999) can be computed from the phase difference of the signals. Tass et al. propose two methods to compute the synchronization index, one based on the Shannon entropy of the distribution of the phase differences and the other based on the conditional probability of the phases. The integer values n and m in this method are determined by trial and error, i.e. the synchronization index is computed for many combinations of n and m and those numbers which yield the largest index are selected for further analysis. Lachaux et al.

define the PLV as PLVt = 1/N

where the sum is taken over the number of trials. The PLV is also computed for surrogate data generated by shuffling the trial indices for one signal, while keeping the other constant, which allows for assessing the statistical significance of the observed synchronization. A comparison has shown that there is no fundamental difference between the synchronization index and PLV and that they are equally well suited for detecting phase synchrony (Le Van Quyen, et al., 2001). If applied in the channel domain, as noted by Lauchaux, et al., 1999, both methods also suffer from the crosstalk effect and can produce false synchrony measures, even in the absence of synchronous sources. Since sensors pick up signals from multiple sources and the phases of these sources do not combine linearly, crosstalk can result in a nonzero-phase difference, thus making it hard to distinguish from genuine phase locking.

Analysis of phase locking has revealed information about cortico-cortical and cortico-muscular coupling of signals (Fell, et al., 2001, Gross, et al., 2001, Gross, et al., 2000, Lauchaux, et al., 1999, Tass, et al., 1998). The limitation of phase locking methods is that they typically only consider interactions within a narrow frequency band. True large scale neural interactions may involve interactions between different frequencies, requiring consideration of a larger bandwidth (Lauchaux, et al., 1999, Tallon-Baudry, et al., 1997). An example showing the PLV between left and right motor cortices is shown in Fig. 3.

Fig. 3. Example of a PLV map computed using the wavelet transform for frequencies between 3 and 30 Hz. Phase locking was computed for the time series of the left and right motor cortices over a period of 400 ms prior to a movement of the right index finger. The time series were estimated from 275 channel MEG data using the minimum norm inverse

Fig. 3. Example of a PLV map computed using the wavelet transform for frequencies between 3 and 30 Hz. Phase locking was computed for the time series of the left and right motor cortices over a period of 400 ms prior to a movement of the right index finger. The time series were estimated from 275 channel MEG data using the minimum norm inverse

Higher Order Spectral Analysis

Higher order spectral analysis (Nikias and Mendel, 1993) can provide information about non-linear coupling between frequencies of different signals or within a single signal and can be seen as an extension of the classical power spectrum. While the power spectrum can be computed from the Fourier transform of the autocorrelation of the signal (and likewise, the cross spectrum from the Fourier transform of the cross correlation of two signals), higher order spectra are computed from the third, fourth or nth order cumulants of the signals. The Bispectrum for example is computed as the two-dimensional Fourier transform of the third-order cumulant of the signal. Let x(k) be a nth order stationary discrete time series, then the nth order cumulants are defined as

Cn (ti,T2, . . .,rn-i) = mn (TiT. .. ,T—i) - m% (TiT. . . ,T—i) where mn (ti,T2, . . .,Tn-l) = E [x (k) X (k + Ti) . . .X (k + Tn-1)] , is the nth moment of x, and m% is the respective moment of an equivalent Gaussian signal with the same mean and autocorrelation as x (Nikias and

Mendel, 1993). The nth order spectrum is thendefined as:

oo oo

Pn (/1,/2,...,/n-1) = e ••• Cn (t1 ,t2,-.,Tn-1)

Vk=1

In a similar fashion, the nth order cross spectra can be computed by using the Fourier transforms of the cross-cumulants. For example, the bispectrum P3 (fi, f2) between three signals is defined as the 2D Fourier transform of the cross-cumulants:

4vz (T!, t2 ) = E [x (k) y (k + n) z (k + r2)j - mfyz (n, r2)

The cross-bicoherence can be computed by normalizing by the power spectral densities of the component signals:

Cvz (fi, f2) = P3 (/1./2) . The bispectrum and cross-

bispectrum, as well as the bicoherence and cross-bicoherence, can detect quadratic phase coupling between signals x, y and z. By replacing z with either x or y, we can use the cross bispectrum to detect coupling between any pair of signals as illustrated in Fig 4. Due to the asymmetry of the cross-bispectrum

Fig. 4. Sample bicoherence map for the same time series described in Fig. 3. The map shows a non-linear interaction of 12 Hz and 20 Hz prior to the finger movement between the left and right motor cortex

for two signals, i.e. C%xy (fi, f2) = Cyxx (fi, f2), this measure can potentially be used to infer the direction of interaction from one signal to another.

Higher order spectra have many desirable properties that make them useful tools for analyzing neuronal interaction. The method is not limited to narrowband signals and the higher order spectra reveal phase information, which makes the method well suited to the analysis of time phase coupling (Jamsek, et al., 2003, Schack, et al., 2002). Because the higher order moments for a signal with a Gaussian distribution vanish, the higher order spectra are insensitive to Gaussian noise. Furthermore, in contrast to the coherence metric, the bispectrum is less sensitive to linear coupling and therefore is an attractive alternative for the detection of nonlinear interactions (Jamsek, et al., 2003). Applications of higher order spectral analysis to EEG data in the channel domain have been presented by Pfurtscheller and Lopes de Silva (1999), where non-linear interaction between 11 Hz and 22Hz components in the post-movement beta event-related synchronization (ERS) was demonstrated, as well as in the analysis of short-term memory (Schack, et al., 2002) and in microelectrode recordings from the visual cortices of cats and monkeys (Schanze and Eckhorn, 1997). Because of its insensitivity to linear relationships and therefore its robustness to linear cross talk effects between signals, the method is also well suited for application in the source domain.

Structural Equation Modeling

While the methods described so far in this chapter can be used to analyze the interactions between two signals, it is also of great interest to perform network analysis involving multiple sources or ensembles of neuronal activity. A simple approximation to the potentially complex interactions between multiple regions is given by the structural equation model (SEM) or path analysis, which assumes a linear relationship between the activity in each region and also with respect to any experimental variables (Astolfi, et al., 2005, Friston, 1994,, Bollen 1989). This assumption can be cast in a simple equation (Mcintosh and Gonzalez-Lima, 1994):

y = B • y+£■ x + z y,z € Rm,x € Rn, B € RnXn, r € Rmxn (4)

The vector y represents the signal at a single time point for each of the m connected components of the network and the vector x represents n independent external components, e.g. variations of the experimental conditions. These can be binary variables that are set to one for the experimental condition and zero for the control condition (Mcintosh and Gonzalez-Lima, 1994). The matrices B and r represent, respectively, the influence of the network on itself and the influence of experimental conditions on the network. In SEM it is assumed that the components of the network do not interact with themselves, therefore the diagonal elements of the matrix B are set to zero. The vector Z models the residual activity not explained by the linear model. The model parameters are estimated from the observed signals using least squares regression (Astolfi et al., 2005). Since typically the number of unknowns in this optimization problem exceeds the number of equations, a priori information about the model has to be provided. This is done by limiting the number of free parameters, i.e. the elements of the matrix B. By setting some elements of B to zero, interactions between selected components are precluded and the number of unknowns is reduced, making the system solvable. Because the method is based on the covariance it will be affected by linear crosstalk leading to potentially erroneous inferences of network connectivity when applied to EEG or MEG data. The SEM approach was originally applied in functional neuroimaging to PET and fMRI studies in which dynamic data were not available. However, when applied to dynamic EEG and MEG data, the assumption of instantaneous linear unidirectional interactions between the network components is rather restrictive. The use of only the zero-lag co-variance between regions ignores the temporal structure of the data, so that shuffling the data in time would have no influence on the SEM parameters (Ramnani, et al., 2004). We conclude this chapter with a brief introduction to two network models that do consider dynamic interactions: the linear MVAR and the nonlinear dynamic causal models.

Multivariate Autoregressive Models (MVAR)

Multivariate autoregressive models can also be used to model interactions between multiple regions. However, unlike SEM, they do not assume instantaneous linear interaction between the regions, but also take the past of each signal into account. One of the key attractions of the MVAR model is its ability to make inferences about the direction of interaction of multiple network components (Brovelli, et al., 2004). The mth order MVAR model is described by the following equation: