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 subjects brain. Furthermore, these distributions of source power can be subjected to statistical analysis.

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 ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/beamformer/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);

• Why is it important that the length of each data piece is the length of a fixed number of oscillatory cycles?

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 segmentedmriF'.

Otherwise, segmentation involves the following steps ¹⁾:

mri = ft_read_mri('Subject01.mri');
cfg = [];
cfg.write        = 'no';
cfg.coordinates  = 'ctf';
cfg.template     = '/PATH/TO/spm2/templates/T1.mnc'; % at the Donders, you don't need to specify this
[segmentedmri] = ft_volumesegment(cfg, mri);

The following commands pertains to a bug not yet fixed. Quite often it happens that the segmented volumes don't align anymore with the anatomy. You can check this by

figure;imagesc(mri.anatomy(:,:,100));
figure;imagesc(segmentedmri.gray(:,:,100));

The following lines serve to flip the segmented brain in order to fit the anatomy:

segmentedmriF = segmentedmri;
segmentedmriF.gray  = flipdim(flipdim(flipdim(segmentedmriF.gray,3),2),1);
segmentedmriF.white = flipdim(flipdim(flipdim(segmentedmriF.white,3),2),1);
segmentedmriF.csf   = flipdim(flipdim(flipdim(segmentedmriF.csf,3),2),1);

Save the segmented MRI:

save segmentedmriF segmentedmriF

If you don't have SPM installed, but nevertheless want to continue with this tutorial, the segmented mri is available from ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/beamformer/segmentedmriF.mat. You can load it from disk and plot it to check that the segmented brain fits and the MRI are aligned:

load segmentedmriF

mri = ft_read_mri('Subject01.mri');
segmentedmriF.transform = mri.transform;
segmentedmriF.anatomy   = mri.anatomy;

figure
cfg = [];
ft_sourceplot(cfg,segmentedmriF); %only mri
figure
cfg.funparameter = 'gray';
ft_sourceplot(cfg,segmentedmriF); %segmented gray matter on top
figure
cfg.funparameter = 'white';
ft_sourceplot(cfg,segmentedmriF); %segmented white matter on top
figure
cfg.funparameter = 'csf';
ft_sourceplot(cfg,segmentedmriF); %segmented csf matter on top

Now prepare the head model from the segmented brain surface:

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

• Why might a single sphere model be inadequate for performing beamformer estimates?

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. Important: sensors MLP31 and MLO12 were removed from the data set. Thus it is essential to remove these sensors as well when calculating the lead fields.

cfg                 = [];
cfg.grad            = freqPre.grad;
cfg.vol             = vol;
cfg.reducerank      = 2;
cfg.channel         = {'MEG','-MLP31', '-MLO12'};
cfg.grid.resolution = 1;
[grid] = ft_prepare_leadfield(cfg);

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:

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 the output:

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 ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/beamformer/source.mat) and the subject's MRI (which is available from ftp://ftp.fcdonders.nl/pub/fieldtrip/tutorial/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 feature of co-variance based spatial filters attributing more power to deep sources.

• Discuss why the source power is overestimated in the center of the brain

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.coordinates   = '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

cfg = [];
cfg.coordinates  = 'ctf';
cfg.nonlinear  = 'no';
sourceDiffIntN = ft_volumenormalise(cfg, sourceDiffInt);

Now the data can be plotted

sourceDiffIntN = rmfield(sourceDiffIntN,'inside'); %due to a bug the 'inside' field has to be removed first

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, sourceDiffIntN);
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 sourceanalysis. This was done on the basis of the smallest eigen value 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.

• Compare figure 3 and 7. It appears that normalizing the power with the baseline activity result in fewer and more focal sources. Why?

• On page 5 the regularization parameter was cfg.lambda = 0. Change it to cfg.lambda = 1e-29 and plot the power estimate with respect to baseline. How does the regularization parameter affect the properties of the spatial filter?

This tutorial has last been tested with version 20100329 of FieldTrip