Coherence

The structure of this tutorial does not yet follow the documentation guidelines.

To study the oscillatory synchrony between two signals, one can compute the coherence. This is computed in the frequency domain by normalizing the magnitude of the summed cross-spectral density between two signals by their respective power. For each frequency bin the coherence value is a number between 0 and 1. The coherence values reflect the consistency of the phase difference between the two signals at a given frequency. In this session we will explore the concept of coherence by investigating a dataset from an experiment in which the subject was required to maintain an isometric contraction of a forearm muscle. The coherence between the MEG signals and the acquired EMG will be estimated. First we will explore the coherence between the EMG signal and all MEG channels. Secondly, we will investigate how the coherence estimate is influenced by the number of trials, and by the degree of spectral smoothing using multitaper spectral analysis. Even though the example in this session covers cortico-muscular coherence, coherence between sensors can be calculated in exactly the same way.

To compute the coherence between the MEG and EMG signals for the example dataset we will perform the following steps:

• Read the data into Matlab using ft_preprocessing
• Compute the power spectra and cross-spectral densities using the function ft_freqanalysis and subsequently compute the coherence using ft_connectivityanalysis
• Visualize the results using ft_singleplotER, ft_multiplotER, and ft_topoplotER
• Subsequently it is possible to localise the neuronal sources coherent with the EMG, using ft_sourceanalysis

The dataset used in this example has been recorded in an experiment in which the subject had to lift her hand and exert a constant force against a lever. The force was monitored by strain gauges on the lever. The subject performed two blocks of 20 trials in which either the left or the right wrist was extended for about 10 seconds. A trial started as soon as the subject managed to get his force output within a specified range from 1 to 2 N. If the force was not kept constant during the course of a trial, the trial was terminated prematurely.

The bipolar EMG signal was recorded from the right extensor carpi radialis longus muscle in the lower arm. MEG signals were recorded with a 151 sensor CTF Omega System (Port Coquitlam, Canada). In addition, the EOG was recorded to later discard trials contaminated by eye movements and blinks. The ongoing MEG and EOG signals were lowpass filtered at 300 Hz, digitized at 1200 Hz and stored for off-line analysis. To measure the head position with respect to the sensors, three coils were placed at anatomical landmarks of the head (nasion, left and right ear canal). While the subjects were seated under the MEG helmet, the positions of the coils were determined before and after the experiment by measuring the magnetic signals produced by currents passed through the coils. Magnetic resonance images (MRIs) were obtained from a 1.5 T Siemens system. During the MRI scan, ear molds containing small containers filled with vitamin E marked the same landmarks. This allows us, together with the anatomical landmarks, to align source estimates of the MEG with the MRI.

We will calculate the coherence between the MEG and the EMG when the subject extended her LEFT wrist, while keeping the right forearm muscle relaxed. The first step is to read the data. In this section we will apply automatic artifact rejection. Preprocessing requires the original MEG dataset, which is available from ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/SubjectCMC.zip.

The epochs of interest have to be defined according to a custom-written function called trialfun_left.m. Note that this function is not part of the FieldTrip toolbox: see appendix 2, or download it from ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/coherence/trialfun_left.m. This function uses the information provided by the triggers which were recorded simultaneously with the data. In this experiment each trigger corresponds with the start or the end of a contraction. The epochs which correspond to a contraction of the left forearm muscle are selected. Subsequently these 10-second pieces are cut into ten 1-second epochs.

% find the interesting epochs of data
cfg = [];
cfg.trialfun                  = 'trialfun_left';
cfg.dataset                   = 'SubjectCMC.ds';
cfg = ft_definetrial(cfg);

% detect EOG artifacts in the MEG data
cfg.continuous                = 'yes';
cfg.artfctdef.eog.padding     = 0;
cfg.artfctdef.eog.bpfilter    = 'no';
cfg.artfctdef.eog.detrend     = 'yes';
cfg.artfctdef.eog.hilbert     = 'no';
cfg.artfctdef.eog.rectify     = 'yes';
cfg.artfctdef.eog.cutoff      = 2.5;
cfg.artfctdef.eog.feedback    = 'no';
cfg = ft_artifact_eog(cfg);

% detect jump artifacts in the MEG data
cfg.artfctdef.jump.feedback   = 'no';
cfg.padding                   = 5;
cfg = ft_artifact_jump(cfg);

% detect muscle artifacts in the MEG data
cfg.artfctdef.muscle.cutoff   = 8;
cfg.artfctdef.muscle.feedback = 'no';
cfg = ft_artifact_muscle(cfg);

% reject the epochs that contain artifacts
cfg.artfctdef.reject          = 'complete';
cfg = ft_rejectartifact(cfg);

% preprocess the MEG data
cfg.blc                       = 'yes';
cfg.dftfilter                 = 'yes';
cfg.channel                   = {'MEG'};
meg = ft_preprocessing(cfg);

Next, read the left and right EMG data. Note that the settings are different for the EMG and MEG data. Most importantly, the EMG data are highpass filtered and rectified. This is a standard procedure when calculating cortico-muscle coherence.

cfg              = [];
cfg.dataset      = meg.cfg.dataset;
cfg.trl          = meg.cfg.trl;
cfg.continuous   = 'yes';
cfg.blc          = 'yes';
cfg.dftfilter    = 'yes';
cfg.channel      = {'EMGlft' 'EMGrgt'};
cfg.hpfilter     = 'yes';
cfg.hpfreq       = 10;
cfg.rectify      = 'yes';
emg = ft_preprocessing(cfg);

Finally, combine the EMG and MEG trials to a common data structure:

data = ft_appenddata([], meg, emg);

Note, that due to the artifact rejection, this procedure is very slow and we typically would want to perform this step only once. Therefore you can save the preprocessed data:

save data data 

To get a feel for the data, plot a trial from a sensor overlying the left motor-cortex (MRC21) and the left and right EMG-signals, by selecting the first trial from the data:

figure
subplot(2,1,1);
plot(data.time{1},data.trial{1}(77,:));
axis tight;
legend(data.label(77));

subplot(2,1,2);
plot(data.time{1},data.trial{1}(152:153,:));
axis tight;
legend(data.label(152:153));

Figure 1; An example of the raw MEG data from sensor MLC21 (upper frame) and the EMG data (lower frame). The signals are from the output of ft_preprocessing and plotted using the matlab plot function. Note that the signal strength of the left EMG is bigger than that of the right EMG.

Explore the MEG and EMG in figure 1, e.g. by zooming in. How are the signals different from one another?

Using ft_freqanalysis, the characteristics in the frequency domain will be computed. This step requires the preprocessed MEG and EMG data (see above or download from ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/coherence/data.mat). Load the data with:

load data 

First, a configuration structure (cfg) must be defined. The FFT-algorithm will be used to compute the fourier representation of each signal. To optimize the estimation, spectral smoothing using ‘multitapers’ will be applied. In this context, the degree of ‘smoothing’ (as defined in the parameter cfg.tapsmofrq) is critical. We will return to this parameter later.

cfg            = [];
cfg.output     = 'powandcsd';
cfg.method     = 'mtmfft';
cfg.foilim     = [5 100];
cfg.tapsmofrq  = 5;
cfg.keeptrials = 'yes';
cfg.channel    = {'MEG' 'EMGlft' 'EMGrgt'};
cfg.channelcmb = {'MEG' 'EMGlft'; 'MEG' 'EMGrgt'};
freq           = ft_freqanalysis(cfg, data);

To calculate the coherence between the EMG and the MEG signals from the powerspectra and cross-spectral densities use the following function:

cfg            = [];
cfg.method     = 'coh';
cfg.channelcmb = {'MEG' 'EMG'};
fd             = ft_connectivityanalysis(cfg, freq);

The field fd.cohspctrm now contains the coherence for all MEG sensors with respect to the EMG signals.

Visualize the coherence between the EMG and all the MEG sensors:

cfg                  = [];
cfg.xparam           = 'freq';
cfg.zparam           = 'cohspctrm';
cfg.xlim             = [5 80];
cfg.cohrefchannel    = 'EMGlft';
cfg.layout           = 'CTF151.lay';
cfg.showlabels       = 'yes';
figure; ft_multiplotER(cfg, fd)

Figure 2; The coherence between the left EMG and all the MEG sensors calculated using ft_freqanalysis and ft_connectivityanalysis. Plotting was done with ft_multiplotER.

Plot the coherence for sensor MRC21 (using the same settings as in ft_multiplotER):

cfg.channel = 'MRC21';
figure; ft_singleplotER(cfg, fd);

Figure 3; The coherence spectrum between the EMG and sensor MRC21.

a) What determines the frequency resolution of the spectrum, as displayed in figure 3? How can it be increased or decreased? Answer the same question for smoothing.

b) Plot a topographical distribution of the coherence in the beta band. The variable cfg.xlim defines the edges of the frequency band.

cfg                  = [];
cfg.xparam           = 'freq';
cfg.zparam           = 'cohspctrm';
cfg.xlim             = [15 20];
cfg.zlim             = [0 0.1];
cfg.cohrefchannel    = 'EMGlft';
cfg.layout           = 'CTF151.lay';
figure; ft_topoplotER(cfg, fd)

Figure 4; A topographic representation of the coherence between the left EMG and the sensors. The plot was created with ft_topoplotER.

a) Explain the pattern of activation in Figure 4.

b) Plot the topographic representation for other frequencies that might be of interest.

Explore the consequence of changing the smoothing in the frequency domain. Do this by recomputing the cortico-muscular coherence between the EMG signal and MEG sensor MRC21 for different degrees of smoothing. Compute the powerspectra and the cross-spectra, and the corresponding coherence using different degrees of smoothing.

a) 2 Hz smoothing (cfg.tapsmofrq = 2 Hz)

cfg            = [];
cfg.output     = 'powandcsd';
cfg.method     = 'mtmfft';
cfg.foilim     = [5 100];
cfg.tapsmofrq  = 2;
cfg.keeptrials = 'yes';
cfg.channel    = {'MEG' 'EMGlft'};
cfg.channelcmb = {'MEG' 'EMGlft'};
freq2          = ft_freqanalysis(cfg,data);

cfg            = [];
cfg.method     = 'coh';
cfg.channelcmb = {'MEG' 'EMG'};
fd2            = ft_connectivityanalysis(cfg,freq2);

Plot the results of the 5 and 2Hz smoothing:

cfg               = [];
cfg.xparam        = 'freq';
cfg.zparam        = 'cohspctrm';
cfg.cohrefchannel = 'EMGlft';
cfg.xlim          = [5 80];
cfg.channel       = 'MRC21';
figure; ft_singleplotER(cfg, fd, fd2);

b) 10 Hz smoothing (e.g. cfg.tapsmofrq = 10 Hz)

cfg            = [];
cfg.output     = 'powandcsd';
cfg.method     = 'mtmfft';
cfg.foilim     = [5 100];
cfg.keeptrials = 'yes';
cfg.channel    = {'MEG' 'EMGlft'};
cfg.channelcmb = {'MEG' 'EMGlft'};
cfg.tapsmofrq = 10;
freq10        = ft_freqanalysis(cfg,data);

cfg            = [];
cfg.method     = 'coh';
cfg.channelcmb = {'MEG' 'EMG'};
fd10          = ft_connectivityanalysis(cfg,freq10);

Plot the results of the 5, 2, and 10 Hz smoothing:

cfg               = [];
cfg.xparam        = 'freq';
cfg.zparam        = 'cohspctrm';
cfg.xlim          = [5 80];
cfg.ylim          = [0 0.2];
cfg.cohrefchannel = 'EMGlft';
cfg.channel       = 'MRC21';
figure;ft_singleplotER(cfg, fd, fd2, fd10);

Which degree of smoothing do you consider optimal in the calculations above?

Another question pertains to how the estimate of coherence is affected by the number of trials. We will compare the cortico-muscular coherence at two MEG sensors for different amount of data.

Create the following configuration, and compute the coherence.

cfg            = [];
cfg.output     = 'powandcsd';
cfg.method     = 'mtmfft';
cfg.foilim     = [5 100];
cfg.tapsmofrq  = 5;
cfg.keeptrials = 'yes';
cfg.channel    = {'MEG' 'EMGlft'};
cfg.channelcmb = {'MEG' 'EMGlft'};
cfg.trials     = 1:50;  
freq50         = ft_freqanalysis(cfg,data);

cfg            = [];
cfg.method     = 'coh';
cfg.channelcmb = {'MEG' 'EGM'};
fd50           = ft_connectivityanalysis(cfg,freq50);

Plot the results:

cfg                  = [];
cfg.xparam           = 'freq';
cfg.zparam           = 'cohspctrm';
cfg.xlim             = [5 100];
cfg.ylim             = [0 0.2];
cfg.cohrefchannel    = 'EMGlft';
cfg.channel          = 'MRC21';
figure; ft_singleplotER(cfg, fd, fd50);

Compare the results with figure 3. Pay special attention to the noise bias.

In order to localise the neuronal sources which are coherent with the EMG, we can apply beamformers to the data. For a more extensive background in beamforming, in particular beamforming with frequency-domain data, please consult the beamformer tutorial. In this example, we are going to use an algorithm, known as DICS, to estimate the activity of the neuronal sources and to subsequently estimate the coherence with the EMG. In order to achieve this, we first need an estimate of the cross-spectral density between all MEG-channel combinations, and between the MEG-channels and the EMG, at a frequency of interest. This requires the preprocessed data, see above, or download from ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/coherence/data.mat. Load with:

load data

Compute the cross-spectral density matrix for 18 Hz:

cfg            = [];
cfg.method     = 'mtmfft';
cfg.output     = 'powandcsd';
cfg.foilim     = [18 18];
cfg.tapsmofrq  = 5;
cfg.keeptrials = 'yes';
cfg.channelcmb = {'MEG' 'MEG';'MEG' 'EMGlft'};
freq           = ft_freqanalysis(cfg, data);

Once we computed this, we can use ft_sourceanalysis using the following configuration.
This step requires the subject's headmodel, which is available from ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/SubjectCMC.zip.

cfg           = [];
cfg.method    = 'dics';
cfg.refchan   = 'EMGlft';
cfg.frequency = 18;
cfg.hdmfile   = 'SubjectCMC.hdm';
cfg.inwardshift     = 1;
cfg.grid.resolution = 1;
source        = ft_sourceanalysis(cfg, freq);

The resulting source-structure is a volumetric reconstruction which is specified in head-coordinates. In order to be able to visualise the result with respect to the subject's MRI, we have to interpolate the functional data to the anatomical MRI.
For this, we need the subject's MRI, which is available from ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/SubjectCMC.zip.

mri = ft_read_mri('SubjectCMC.mri');

Next, we can proceed with the interpolation.

cfg            = [];
cfg.parameter  = 'coh';
cfg.downsample = 2;
interp         = ft_sourceinterpolate(cfg, source, mri);

There are various ways to visualise the volumetric interpolated data. The most straightforward way is using ft_sourceplot.

cfg              = [];
cfg.method       = 'ortho';
cfg.interactive  = 'yes';
cfg.funparameter = 'coh';
figure; ft_sourceplot(cfg, interp);

Figure 5; The neuronal source showing maximum coherence with the left EMG at 18 Hz. The plot was created with ft_sourceplot.

The trial function used to extract the trials:

function trl = trialfun_left(cfg)

% read in the triggers and create a trial-matrix
% consisting of 1-second data segments, in which 
% left ECR-muscle is active.

event = ft_read_event(cfg.dataset);
trig  = [event(find(strcmp('backpanel trigger', {event.type}))).value];
indx  = [event(find(strcmp('backpanel trigger', {event.type}))).sample];

%left-condition
sel = [find(trig==1028):find(trig==1029)];

trig = trig(sel);
indx = indx(sel);
 
trl = [];
for j = 1:length(trig)-1
  trg1 = trig(j);
  trg2 = trig(j+1);
  if trg1<=100 & trg2==2080,
    trlok      = [[indx(j)+1:1200:indx(j+1)-1200]' [indx(j)+1200:1200:indx(j+1)]'];
    trlok(:,3) = [0:-1200:-1200*(size(trlok,1)-1)]';
    trl        = [trl; trlok];
  end
end

This tutorial was last tested with version 20100329 of FieldTrip.