# 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 will learn about applying beamformer techniques in the frequency domain. You will learn how to compute appropriate time-frequency windows, an appropriate head model and lead field matrix, and various options for contrasting the effect of interest against some control/baseline. Finally, you will be shown several options for plotting the results overlaid on a structural MRI.

It is expected that you understand the previous steps of preprocessing and filtering the sensor data. Some understanding of the options for computing the head model and forward lead field is also useful.

This tutorial will not cover the time-domain option for LCMV/SAM beamformers (described in Background), nor for beamformers applied to evoked/averaged data (although see an example of how to calculate virtual sensors using LCMV for an example of this).

This tutorial contains hands-on material that we use for the MEG/EEG toolkit course and is complemented by this lecture.

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 inverse filter is based on minimizing the source power (or variance) at a given location, subject to 'unit-gain constraint'. This latter part means that, if a source had power of amplitude 1 and was projected to the sensors by the lead field, the inverse filter applied to the sensors should then reconstruct power of amplitude 1 at that location. Beam forming assumes that sources in different parts of the brain are not temporally correlated.

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 (DICS) and the estimates are calculated in the frequency domain (Gross ET al. 2001). Other beam-former 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:

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 1).

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.

Ft_definetrial and ft_preprocessing require the original MEG dataset, which is available from ftp://ftp.fieldtriptoolbox.org/pub/fieldtrip/tutorial/Subject01.zip.

```% find the interesting segments of data
cfg = [];                                           % empty configuration
cfg.dataset                 = 'Subject01.ds';       % name of CTF dataset
cfg.trialdef.eventtype      = 'backpanel trigger';
cfg.trialdef.prestim        = 1;
cfg.trialdef.poststim       = 2;
cfg.trialdef.eventvalue     = 3;                    % event value of FIC
cfg = ft_definetrial(cfg);

% remove the trials that have artifacts from the trl
cfg.trl([15, 36, 39, 42, 43, 49, 50, 81, 82, 84],:) = [];

% preprocess the data
cfg.channel   = {'MEG', '-MLP31', '-MLO12'};        % read all MEG channels except MLP31 and MLO12
cfg.demean    = 'yes';                              % do baseline correction with the complete trial

dataFIC = ft_preprocessing(cfg);```

These data have been cleaned from artifacts by removing several trials and two sensors; see the visual artifact rejection tutorial.

Subsequently you can save the data to disk.

`save dataFIC dataFIC`
2009/03/11 13:41

### Time windows of interest

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.

1. Suitable control windows are, for example:
1. Activity contrasted with baseline (example shown here using dataPre)
2. Activity of condition 1 contrasted with condition 2 (example not shown)
2. However, if no other suitable data condition or baseline time-window exists, then
1. Activity contrasted with estimated noise (example shown below)
2. Use normalized leadfields (mentioned in 'the forward model and lead field matrix' section and Exercise 4 below)

The null hypothesis for both options within (1) 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.

### Exercise 1: data length

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. The DICS 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. It is given as output when cfg.output = 'powandcsd'. 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 a given 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 1):

```mri = ft_read_mri('Subject01.mri');
cfg = [];
cfg.write      = 'no';
[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 = [];
cfg.method = 'singleshell';
If you want to do a beamformer source reconstruction on EEG data, you have to pay special attention to the EEG referencing. The forward model will be made with an common average reference [*], i.e. the mean value over all electrodes is zero. Consequently, this also has to be true in your data.

Prior to doing the spectral decomposition with ft_freqanalysis you have to ensure with ft_preprocessing that all channels are re-referenced to the common average reference.

Furthermore, after selecting the channels you want to use in the sourcereconstruction (excluding the bad channels) and after re-referencing them, you should not make sub-selections of channels any more and throw out channels, because that would cause the data not be average referenced any more.

[*] except in some rare cases, like with bipolar iEEG electrode montages

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.

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.reducerank      = 2;
cfg.channel         = {'MEG','-MLP31', '-MLO12'};
cfg.grid.resolution = 1;   % use a 3-D grid with a 1 cm resolution
cfg.grid.unit       = 'cm';

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 help 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. We first show the example as if you only had dataPost (and no dataPre against which to make a contrast).

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

sourcePost_nocon = ft_sourceanalysis(cfg, freqPost);```

The purpose of cfg.dics.projectnoise will become more clear in the section on Neural Activity Index. The purpose of lambda is discussed in Exercise 6.

Save the output:

`save sourcePost_nocon sourcePost_nocon`

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 (sourcePost_nocon.mat)) and the subject's MRI (also is available from the ftp server (Subject01.zip)).

`load sourcePost_nocon`

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. Using the ft_volumereslice function before doing the interpolation ensures that the MRI is well behaved, because the reslicing causes the voxel axes to be aligned with the head coordinate axes:

```mri = ft_read_mri('Subject01.mri');
mri = ft_volumereslice([], mri);

cfg            = [];
cfg.downsample = 2;
cfg.parameter = 'avg.pow';
sourcePostInt_nocon  = ft_sourceinterpolate(cfg, sourcePost_nocon , mri);```

Plot the interpolated data:

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

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. There are several ways of circumventing the noise bias towards the center of the head which we will show below.

### Exercise 3: center of head bias

Discuss why the source power is overestimated in the center of the brain. Hint 1: what are the leadfield values in the center of the head? Why? Hint 2: Remember the 'unit-gain constraint' of beamformer spatial filters.

### Neural Activity Index

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), in order to remove the center of the head bias shown above. 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, by setting cfg.dics.projectnoise='yes' (default is 'no'). This noise estimate was computed on the basis of the smallest eigenvalue of the cross-spectral density matrix. To calculate the NAI do the following:

```sourceNAI = sourcePost_nocon;
sourceNAI.avg.pow = sourcePost_nocon.avg.pow ./ sourcePost_nocon.avg.noise;

cfg = [];
cfg.downsample = 2;
cfg.parameter = 'avg.pow';
sourceNAIInt = ft_sourceinterpolate(cfg, sourceNAI , mri);```

Plot it:

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

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

### Exercise 4: lead field normalization

Another option, besides contrasting to the noise estimate, is to normalize the lead field when you compute it (cfg.normalize='yes' in the call to ft_prepare_leadfield). Recompute the lead field and source estimate this way and plot the result.

One approach is to compare the post- and pre-stimulus interval, which we describe now here. Another approach is to contrast the same window of interest (relative in time to some stimulus or response) between two or more conditions, which we do not show here, but the calls to FieldTrip functions are conceptually the same.

Importantly, if you later want to compare the two conditions statistically, you have to compute the sources based on an inverse filter computed from both conditions, so called 'common filters', and then apply this filter separately to each condition to obtain the source power estimate in each condition separately.

First we compute a single data structure with both conditions, and compute the frequency domain CSD.

`dataAll = ft_appenddata([], dataPre, dataPost);`
```cfg = [];
cfg.method    = 'mtmfft';
cfg.output    = 'powandcsd';
cfg.tapsmofrq = 4;
cfg.foilim    = [18 18];
freqAll = ft_freqanalysis(cfg, dataAll);```

Then we compute the inverse filter based on both conditions. Note the use of cfg.keepfilter so that the output saves this computed filter.

```cfg              = [];
cfg.method       = 'dics';
cfg.frequency    = 18;
cfg.grid         = grid;
cfg.dics.projectnoise = 'yes';
cfg.dics.lambda       = '5%';
cfg.dics.keepfilter   = 'yes';
cfg.dics.realfilter   = 'yes';
sourceAll = ft_sourceanalysis(cfg, freqAll);```

By placing this pre-computed filter inside cfg.grid.filter, it can now be applied to each condition separately.

```cfg.grid.filter = sourceAll.avg.filter;
sourcePre_con  = ft_sourceanalysis(cfg, freqPre );
sourcePost_con = ft_sourceanalysis(cfg, freqPost);```
```save sourcePre_con sourcePre_con
save sourcePost_con sourcePost_con```

Now we can compute the contrast of (post-pre)/pre. 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_con;
sourceDiff.avg.pow = (sourcePost_con.avg.pow - sourcePre_con.avg.pow) ./ sourcePre_con.avg.pow;```

Load and reslice the MRI if not done already in previous step:

```mri = ft_read_mri('Subject01.mri');
mri = ft_volumereslice([], mri);```

Then interpolate the source to the MRI:

```cfg            = [];
cfg.downsample = 2;
cfg.parameter  = 'avg.pow';
sourceDiffInt  = ft_sourceinterpolate(cfg, sourceDiff , mri);```

Now plot the power ratios:

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

Figure 4; sourceplot with method “slice”

### Exercise 5: comparing normalizations

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

### Exercise 6: regularization

The regularization parameter was cfg.dics.lambda = '5%'. Change it to 0 or to '10%' and plot the power estimate with respect to baseline. How does the regularization parameter affect the properties of the spatial filter?

To plot an 'orthogonal cut':

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

Figure 5; 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.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.funcolorlim   = [0.0 1.2];
cfg.opacitylim    = [0.0 1.2];
cfg.opacitymap    = 'rampup';
ft_sourceplot(cfg, sourceDiffIntNorm);```

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

You can also project the power onto a surface using ft_sourceplot. 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.nonlinear     = 'no';
sourceDiffIntNorm = ft_volumenormalise(cfg, sourceDiffInt);```

Now the data can be plotted

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

Figure 7; sourceplot with method “surface”

Beamforming source analysis in the frequency domain with DICS has been demonstrated. An example of how to compute a head model (single shell) and forward lead field was shown. Various options for contrasting the time-frequency window of interest against a control was shown, including against 'noise' (NAI: minimum eigenvalue of the CSD) and against the pre-stimulus window using a 'common filter'. Options at each stage and their influence on the results were discussed, such as lead field normalization and CSD matrix regularization. Finally, options for plotting on slices, orthogonal views, or on the surface were shown.

Details on head models can be found here or here. Computing event-related fields with MNE or LCMV might be of interest. More information on common filters can be found here. If you are doing a group study where you want the grid points to be the same over all subjects, see here. See here for source statistics.

FAQs:

 Can I do combined EEG and MEG source reconstruction? 2009/11/09 15:25 Can I restrict the source reconstruction to the grey matter? 2009/02/17 15:18 Robert Oostenveld How are the different head and MRI coordinate systems defined? 2009/02/17 15:16 Robert Oostenveld How are the Left and Right Pre-Auricular (LPA and RPA) points defined? 2012/08/17 10:23 Robert Oostenveld How can I check whether the grid that I have is aligned to the segmented volume and to the sensor gradiometer? 2012/03/02 15:55 How can I determine the anatomical label of a source? 2010/04/25 07:50 Robert Oostenveld How can I fine-tune my BEM volume conduction model? 2013/12/05 11:01 Robert Oostenveld How can I use OpenMEEG for forward modeling? 2011/09/07 09:57 How can I visualize the different geometrical objects that are needed for forward and inverse computations? 2011/11/10 11:01 Jan-Mathijs Schoffelen how to coregister an anatomical mri with the gradiometer or electrode positions 2010/10/06 16:30 Robert Oostenveld Is it important to have accurate measurements of electrode locations for EEG source reconstruction? 2009/02/17 15:18 Robert Oostenveld What is the conductivity of the brain, CSF, skull and skin tissue? 2010/01/12 10:44 What kind of volume conduction models of the head (head models) are implemented? 2009/02/17 15:15 Robert Oostenveld Where can I find the dipoli command-line executable? 2009/02/26 22:06 Where is the anterior commissure? 2016/07/13 10:51 Robert Oostenveld Why is there a rim around the brain for which the source reconstruction is not computed? 2009/02/17 15:16 Robert Oostenveld Why should I use an average reference for EEG source reconstruction? 2016/11/30 10:54 Robert Oostenveld

Example scripts:

 Align EEG electrode positions to BEM headmodel 2009/03/03 21:52 Check the quality of the anatomical coregistration 2016/11/09 12:04 Robert Oostenveld Combined EEG and MEG source reconstruction 2012/08/09 20:01 Lilla Magyari Common filters in beamforming 2014/09/25 19:55 Compute forward simulated data and apply a beamformer scan 2009/03/03 21:52 Compute forward simulated data and apply a dipole fit 2009/03/03 21:52 Create MNI-aligned grids in individual head-space 2009/07/31 11:49 Doing source-analysis with the created headmodel 2009/07/09 20:11 Example use of the compute_leadfield function 2010/07/21 11:04 Robert Oostenveld Fit a dipole to the tactile ERF after mechanical stimulation 2009/05/20 14:57 How to import data from MNE-Python and FreeSurfer 2018/05/30 15:57 Robert Oostenveld Localizing the sources underlying the difference in event related fields 2015/11/29 10:55 Robert Oostenveld Make leadfields using different headmodels 2009/05/20 14:57 Plotting the result of source reconstructing on a cortical mesh 2010/04/22 09:01 Source statistics 2009/03/05 16:56 Saskia Haegens Testing BEM created EEG lead fields 2009/10/29 16:22 Use your own forward leadfield model in an inverse beamformer computation 2009/03/03 21:58

This tutorial was last tested by Robert with version 20150609 of FieldTrip using MATLAB R2014b on a 64-bit OS X computer.

1)
ft_volumesegment makes use of SPM. The necessary SPM-files are located in fieldtripXXX/external/spm8