Localizing oscillatory sources using beamformer techniques

In this tutorial we will continue working on the dataset described in the preprocessing tutorials. Below we will repeat code to select the trials and preprocess the data as described in the first tutorials (trigger based trial selection, visual artifact rejection).

In this tutorial you can find information about applying beamformer techniques in the frequency domain.

In the Time-Frequency Analysis tutorial we identified strong oscillations in the beta band in a language paradigm. The goal of this section is to identify the sources responsible for producing this oscillatory activity. We will apply a beamformer technique. This is a spatially adaptive filter, allowing us to estimate the amount of activity at any given location in the brain. The brain is divided in a regular three dimensional grid and the source strength for each grid point is computed. The method applied in this example is termed Dynamical Imaging of Coherent Sources and the estimates are calculated in the frequency domain (Gross et al. 2001). Other beamformer methods rely on sources estimates calculated in the time domain, e.g. the Linearly Constrained Minimum Variance (LCMV) and Synthetic Aperture Magnetometry (SAM) methods (van Veen et al., 1997; Robinson and Cheyne, 1997). These methods produce a 3D spatial distribution of the power of the neuronal sources. This distribution is then overlaid on a structural image of the subject's brain. Furthermore, these distributions of source power can be subjected to statistical analysis. It is always ideal to contrast the activity of interest against some control/baseline activity. Options for this will be discussed below, but it is best to keep this in mind when designing your experiment from the start, rather than struggle to find a suitable control/baseline after data collection.

To localize the oscillatory sources for the example dataset we will perform the following steps:

• Read the data into Matlab using ft_definetrial and ft_preprocessing
• Compute the cross-spectral density matrix using the function ft_freqanalysis
• Construct a forward model and lead field matrix using ft_volumesegment, ft_prepare_singleshell and ft_prepare_leadfield
• Compute a spatial filter and estimate the power of the sources using ft_sourceanalysis
• Visualize the results, by first interpolating the sources to the anatomical MRI using ft_sourceinterpolate and plotting this with ft_sourceplot.

The aim is to identify the sources of oscillatory activity in the beta band. From the section time-frequency analysis we have identified 18 Hz as the center frequency for which the power estimates should be calculated. We seek to compare the activation in the post-stimulus to the activation in the pre-stimulus interval. We first use ft_preprocessing and ft_redefinetrial to extract relevant data. It is important that the length of each data piece is the length of a fixed number of oscillatory cycles. Here 9 cycles are used resulting in a 9/18 Hz = 0.5 s time window. Thus, the post-stimulus time-window range between 0.8 to 1.3 s and the pre-stimulus interval between -0.5 to 0.0 s (see Figure 2).

Figure 1; The time-frequency presentation used to determine the time- and frequency-windows prior to beamforming. The squares indicate the selected time-frequency tiles for the pre- and post-response.

Now we select the time windows of interest: the pre- and post stimulus windows. This requires the preprocessed data (see above), which is available from the FieldTrip ftp server (dataFIC.mat). Load the data with the following command:

load dataFIC

Now 'cut' out the pre- and post-stimulus time windows:

cfg = [];                                           
cfg.toilim = [-0.5 0];                       
dataPre = ft_redefinetrial(cfg, dataFIC);
   
cfg.toilim = [0.8 1.3];                       
dataPost = ft_redefinetrial(cfg, dataFIC);

As mentioned in the Background, it is ideal to contrast the activity of interest against some control. Several options to do this are:

1. Activity contrasted with baseline (as shown above with dataPre and below in Example 3)
2. Activity of condition 1 contrasted with condition 2 (example not shown)
3. If no other suitable data condition or baseline time-window exists, then
   1. Activity contrasted with estimated noise (shown below in Example 3)
   2. Use normalized leadfields (mentioned in 'the forward model and lead field matrix' section below)

The null hypothesis for options (1) and (2) is that the data in both conditions are the same, and thus the best spatial filter is the one that is computed using both data conditions together (also known as 'common filters'). This common filter is then applied separately to each condition. The rest of this tutorial proceeds using separate filters per condition, but please keep this common filter in mind.

The beamformer technique is based on an adaptive spatial filter. This spatial filter is derived from the frequency counterpart of the covariance matrix: the cross-spectral density matrix. This matrix contains the cross-spectral densities for all sensor combinations and is computed from the Fourier transformed data of the single trials. The frequency of interest is 18 Hz and the smoothing window is +/-4 Hz:

cfg = [];
cfg.method    = 'mtmfft';
cfg.output    = 'powandcsd';
cfg.tapsmofrq = 4;
cfg.foilim    = [18 18];
freqPre = ft_freqanalysis(cfg, dataPre);

cfg = [];
cfg.method    = 'mtmfft';
cfg.output    = 'powandcsd';
cfg.tapsmofrq = 4;
cfg.foilim    = [18 18];
freqPost = ft_freqanalysis(cfg, dataPost);

The first step in the procedure is to construct a forward model. The forward model allows us to calculate an estimate of the field measured by the MEG sensors for given a current distribution. In MEG analysis a forward model is typically constructed for each subject. There are many types of forward models which to various degrees take the individual anatomy into account. We will here use a semi-realistic head model developed by Nolte (2003). It is based on a correction of the lead field for a spherical volume conductor by a superposition of basis functions, gradients of harmonic functions constructed from spherical harmonics.

The first step in constructing the forward model is to find the brain surface from the subjects MRI. This procedure is termed segmentation. Note that segmentation is quite time consuming. If you have access to the preprocessed file you can skip ahead to 'load segmentedmri'.

Otherwise, segmentation involves the following steps ¹⁾:

mri = ft_read_mri('Subject01.mri');
cfg = [];
cfg.write      = 'no';
cfg.coordsys   = 'ctf';
[segmentedmri] = ft_volumesegment(cfg, mri);

Alternatively, you can load the segmented MRI available from the FieldTrip ftp server (segmentedmri.mat):

load segmentedmri

Now prepare the head model from the segmented brain surface:

cfg = [];
vol = ft_prepare_singleshell(cfg, segmentedmri);

The next step is to discretize the brain volume into a grid. For each grid point the lead field matrix is calculated. It is calculated with respect to a grid with a 1 cm resolution.

cfg                 = [];
cfg.grad            = freqPre.grad;
cfg.vol             = vol;
cfg.reducerank      = 2;
cfg.channel         = {'MEG','-MLP31', '-MLO12'};
cfg.grid.resolution = 1;   % use a 3-D grid with a 1 cm resolution
[grid] = ft_prepare_leadfield(cfg);

As mentioned earlier on, if you are not contrasting the activity of interest against another condition or baseline time-window, then you may choose to normalize the lead field (cfg.normalize='yes'), which will control against the power bias towards the center of the head.

Using the cross-spectral density and the lead field matrices a spatial filter is calculated for each grid point. By applying the filter to the Fourier transformed data we can then estimate the power for the pre- and post-stimulus activity. This results in a power estimate for each grid point. The noise projected to each grid point can also be computed, and is discussed later in the tutorial (below Fig. 6). The purpose of lambda is discussed in Exercise 5.

cfg              = []; 
cfg.frequency    = 18;  
cfg.method       = 'dics';
cfg.projectnoise = 'yes';
cfg.grid         = grid; 
cfg.vol          = vol;
cfg.lambda       = 0;

sourcePre  = ft_sourceanalysis(cfg, freqPre );
sourcePost = ft_sourceanalysis(cfg, freqPost);

save source sourcePre sourcePost

The beamformer procedure estimates the power in the beta frequency band at each grid point in the brain volume. The grid of estimated power values can be plotted superimposed on the anatomical MRI. This requires the output of ft_sourceanalysis (see above or download from the FieldTrip ftp server (source.mat)) and the subject's MRI (also is available from the ftp server (Subject01.zip)).

load source

The function ft_sourceinterpolate aligns the measure of power increase with the structural MRI of the subject. The alignment is done according to the anatomical landmarks (nasion, left and right ear canal) that were both determined in the MEG measurement and in the MRI scan:

mri = ft_read_mri('Subject01.mri');

cfg            = [];
cfg.downsample = 2;
sourcePostInt  = ft_sourceinterpolate(cfg, sourcePost , mri);

Plot the interpolated data:

cfg              = [];
cfg.method       = 'slice';
cfg.funparameter = 'avg.pow';
figure
ft_sourceplot(cfg,sourcePostInt);

Figure 2; The power estimates of the post-stimulus activity only at ~18 Hz. Note the strong noise bias toward the center of the head. The image was done using ft_sourceinterpolate and ft_sourceplot.

Notice that the power is strongest in the center of the brain. This is due to the fact that leadfields in the centre of the head have very small values (why is that?). Due to the 'unit gain constraint' of the beamformer spatial filters, the corresponding spatial filter will have very large values and therefore the beamformer will attribute more power to deep sources.

There are several ways of circumventing the noise bias towards the center of the head.

One approach is to compare the post- and pre-stimulus interval. We will here calculate the relative change in power. In this operation we assume that the noise bias is the same for the pre- and post-stimulus interval and it will thus be removed.

sourceDiff = sourcePost;
sourceDiff.avg.pow = (sourcePost.avg.pow - sourcePre.avg.pow) ./ sourcePre.avg.pow;
  
cfg            = [];
cfg.downsample = 2;
sourceDiffInt  = ft_sourceinterpolate(cfg, sourceDiff , mri);

Now plot the power ratios:

cfg = [];
cfg.method        = 'slice';
cfg.funparameter  = 'avg.pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0.0 1.2];
cfg.opacitylim    = [0.0 1.2]; 
cfg.opacitymap    = 'rampup';  
figure
ft_sourceplot(cfg, sourceDiffInt);

Figure 3; sourceplot with method “slice”

To plot an 'orthogonal cut':

cfg = [];
cfg.method        = 'ortho';
cfg.funparameter  = 'avg.pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0.0 1.2];
cfg.opacitylim    = [0.0 1.2]; 
cfg.opacitymap    = 'rampup';  
figure
ft_sourceplot(cfg, sourceDiffInt);

Figure 4; sourceplot with method “ortho”

The FieldTrip function ft_volumenormalise normalises anatomical and functional volume data to a template anatomical MRI. Spatially aligning the source structures for multiple subjects allows you to compute grandaverages and over subjects statistics. Here we will illustrate the use of volume normalisation for one subject.

cfg = [];
cfg.coordsys      = 'ctf';
cfg.nonlinear     = 'no';
sourceDiffIntNorm = ft_volumenormalise(cfg, sourceDiffInt);

When plotting the orthogonal view it is possible to enter interactive mode by specifying cfg.interactive='yes'. This allows you to 'browse' through the brain volume by specifying the location of cut with a mouse click. To exit the interactive mode press 'q'.

cfg = [];
cfg.method        = 'ortho';
cfg.interactive   = 'yes';
cfg.funparameter  = 'avg.pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0.0 1.2];
cfg.opacitylim    = [0.0 1.2]; 
cfg.opacitymap    = 'rampup';  
figure
ft_sourceplot(cfg, sourceDiffIntNorm);

Figure 5; sourceplot with method “ortho” after volume normalisation

You can also project the power onto a surface. FieldTrip has several surface .mat files available. The surface files are in MNI coordinates, so therefore the volume has to be normalized to match those coordinates. This can be done with the FieldTrip function ft_volumenormalise (see above, as well).

cfg = [];
cfg.coordsys      = 'ctf';
cfg.nonlinear     = 'no';
sourceDiffIntNorm = ft_volumenormalise(cfg, sourceDiffInt);

Now the data can be plotted

cfg = [];
cfg.method         = 'surface';
cfg.funparameter   = 'avg.pow';
cfg.maskparameter  = cfg.funparameter;
cfg.funcolorlim    = [0.0 1.2];
cfg.funcolormap    = 'jet';
cfg.opacitylim     = [0.0 1.2]; 
cfg.opacitymap     = 'rampup';  
cfg.projmethod     = 'nearest'; 
cfg.surffile       = 'surface_l4_both.mat';
cfg.surfdownsample = 10; 
figure
ft_sourceplot(cfg, sourceDiffIntNorm);
view ([90 0])

Figure 6; sourceplot with method “surface”

If it is not possible to compare two conditions (e.g. A versus B or post versus pre) one can apply the neural activity index (NAI). The NAI is the power normalized with an estimate of the spatially inhomogeneous noise. An estimate of the noise has been done by ft_sourceanalysis. This was done on the basis of the smallest eigenvalue of the cross-spectral density matrix. To calculate the NAI do the following:

sourceNAI = sourcePost;
sourceNAI.avg.pow = sourcePost.avg.pow ./ sourcePost.avg.noise;
 
cfg = [];
cfg.downsample = 2;
sourceNAIInt = ft_sourceinterpolate(cfg, sourceNAI , mri);

Plot it:

cfg = [];
cfg.method        = 'slice';
cfg.funparameter  = 'avg.pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [4.0 6.2];
cfg.opacitylim    = [4.0 6.2]; 
cfg.opacitymap    = 'rampup';  
figure
ft_sourceplot(cfg, sourceNAIInt);

Figure 7; The neural activity index (NAI) plotted for the post-stimulus time window normalized with respect to the noise estimate.

This tutorial has last been tested by Jan-Mathijs with version 20110409 of FieldTrip, using Matlab 2009b on a 64-bit Linux platform.