diff --git a/workshop/neuroimaging2-2526.md b/workshop/neuroimaging2-2526.md new file mode 100644 index 000000000..4ce69a474 --- /dev/null +++ b/workshop/neuroimaging2-2526.md @@ -0,0 +1,38 @@ +--- +title: Neuroimaging II - Electrophysiological Methods +tags: [neuroimaging2-2526] +--- + +This course is about the analysis of EEG- and MEG-data. The objective of the course is to introduce the student to the most important analysis methods for this type of data. + +{% include markup/yellow %} +This course is given every year as part of the [Cognitive Neuroscience](https://www.ru.nl/en/education/masters/cognitive-neuroscience-research) research master. + +The version that you are looking at is for the academic year 2025/2026. +{% include markup/end %} + +## Assignments + +This course consists of three parts. Note that some of these parts the assignments are not found here but on [Brightspace](https://brightspace.ru.nl/d2l/home/576894). + +### Signal processing (SPED) + +- [SPED1](/workshop/neuroimaging2-2526/sped1) +- [SPED2](/workshop/neuroimaging2-2526/sped2) +- [SPED3](/workshop/neuroimaging2-2526/sped3) +- [SPED4](/workshop/neuroimaging2-2526/sped4) + +### Source reconstruction (SR) + +- [SR1](/workshop/neuroimaging2-2526/sr1) +- [SR2](/workshop/neuroimaging2-2526/sr2) +- [SR3](/workshop/neuroimaging2-2526/sr3) +- [SR4](/workshop/neuroimaging2-2526/sr4) +- [SR5](/workshop/neuroimaging2-2526/sr5) + +### Statistical testing (SED) + +- [SED1](/workshop/neuroimaging2-2526/sed1) +- [SED2](/workshop/neuroimaging2-2526/sed2) +- [SED3](/workshop/neuroimaging2-2526/sed3) +- [SED4](/workshop/neuroimaging2-2526/sed4) \ No newline at end of file diff --git a/workshop/neuroimaging2-2526/sed1.md b/workshop/neuroimaging2-2526/sed1.md new file mode 100644 index 000000000..5889593d1 --- /dev/null +++ b/workshop/neuroimaging2-2526/sed1.md @@ -0,0 +1,6 @@ +--- +title: SED1 - Statistical testing of electrophysiological data +tags: [neuroimaging2-2425] +--- + +Please look at [Brightspace](https://brightspace.ru.nl/d2l/home/502448) for the assignment. diff --git a/workshop/neuroimaging2-2526/sed2.md b/workshop/neuroimaging2-2526/sed2.md new file mode 100644 index 000000000..f767cb63f --- /dev/null +++ b/workshop/neuroimaging2-2526/sed2.md @@ -0,0 +1,6 @@ +--- +title: SED2 - Statistical testing of electrophysiological data +tags: [neuroimaging2-2425] +--- + +Please look at [Brightspace](https://brightspace.ru.nl/d2l/home/502448) for the assignment. diff --git a/workshop/neuroimaging2-2526/sed3.md b/workshop/neuroimaging2-2526/sed3.md new file mode 100644 index 000000000..fbf6aeb84 --- /dev/null +++ b/workshop/neuroimaging2-2526/sed3.md @@ -0,0 +1,6 @@ +--- +title: SED3 - Statistical testing of electrophysiological data +tags: [neuroimaging2-2425] +--- + +Please look at [Brightspace](https://brightspace.ru.nl/d2l/home/502448) for the assignment. diff --git a/workshop/neuroimaging2-2526/sed4.md b/workshop/neuroimaging2-2526/sed4.md new file mode 100644 index 000000000..44ab00452 --- /dev/null +++ b/workshop/neuroimaging2-2526/sed4.md @@ -0,0 +1,6 @@ +--- +title: SED4 - Statistical testing of electrophysiological data +tags: [neuroimaging2-2425] +--- + +Please look at [Brightspace](https://brightspace.ru.nl/d2l/home/502448) for the assignment. diff --git a/workshop/neuroimaging2-2526/sped1.md b/workshop/neuroimaging2-2526/sped1.md new file mode 100644 index 000000000..ca80e66e5 --- /dev/null +++ b/workshop/neuroimaging2-2526/sped1.md @@ -0,0 +1,6 @@ +--- +title: SPED1 - Signal processing of electrophysiological data +tags: [neuroimaging2-2425] +--- + +Please look at [Brightspace](https://brightspace.ru.nl/d2l/home/502448) for the assignment. diff --git a/workshop/neuroimaging2-2526/sped2.md b/workshop/neuroimaging2-2526/sped2.md new file mode 100644 index 000000000..a5e7ca0a2 --- /dev/null +++ b/workshop/neuroimaging2-2526/sped2.md @@ -0,0 +1,6 @@ +--- +title: SPED2 - Signal processing of electrophysiological data +tags: [neuroimaging2-2425] +--- + +Please look at [Brightspace](https://brightspace.ru.nl/d2l/home/502448) for the assignment. diff --git a/workshop/neuroimaging2-2526/sped3.md b/workshop/neuroimaging2-2526/sped3.md new file mode 100644 index 000000000..5f46def2e --- /dev/null +++ b/workshop/neuroimaging2-2526/sped3.md @@ -0,0 +1,6 @@ +--- +title: SPED3 - Signal processing of electrophysiological data +tags: [neuroimaging2-2425] +--- + +Please look at [Brightspace](https://brightspace.ru.nl/d2l/home/502448) for the assignment. diff --git a/workshop/neuroimaging2-2526/sped4.md b/workshop/neuroimaging2-2526/sped4.md new file mode 100644 index 000000000..abebf7dc1 --- /dev/null +++ b/workshop/neuroimaging2-2526/sped4.md @@ -0,0 +1,352 @@ +--- +title: SPED4 - Time-frequency analysis in practice using FieldTrip +tags: [neuroimaging2-2526] +--- + +## Introduction + +This is an adapted version of the [general FieldTrip tutorial on time-frequency analysis](/tutorial/sensor/timefrequencyanalysis), made specifically for the module _Signal Processing for Electrophysiological Data (SPED)_ in the course _Neuroimaging 2: Electrophysiological Methods_, in the _Cognitive Neuroscience Masters (CNS)_ program at Radboud University. This module is taught by Eelke Spaak. + +Some of the concepts convered here should by now be familiar to you, while some other concepts will be new. The main purpose of this week's tutorial is to show you how spectral and time-frequency data analysis is done when using a typical, widely used, analysis toolbox. This is in contrast to previous weeks, in which the focus was much more on implementing the basic methods yourself and understanding their intricacies. Specifically, this tutorial will cover the time-frequency analysis of a single subject's MEG data using a Hanning window, multitapers and wavelets. The tutorial also shows how to visualize the results. + +### Handing in the hands-on: practicalities + +The results of this hands-on session should be handed in via Brightspace. Please hand in two files: a Matlab script (\*.m) that contains the code you wrote (and/or copied). In that file, please separate all code using code cells (remember `%% code cell header` etc. to split code cells). Each code cell should correspond to one analysis step and/or exercise (if the exercise involves writing/adapting code). In addition to the Matlab script, please copy (figure menu Edit > Copy Figure) the relevant output plots to a separate document (paste in e.g. Word), and write the text-based answers to exercise questions in that document as well. This file should also be handed in (as PDF or Word). + +## Details on the dataset + +The MEG data set used here is from a language study on semantically congruent and incongruent sentences that is described in detail in Wang et al. (2012). Three types of sentences were used in the experiment. In the fully congruent condition (FC) the sentences ended with a high-cloze probability word, e.g., _De klimmers bereikten eindelijk de top van de berg_ (_The climbers finally reached the top of the mountain_) In the fully incongruent condition (FIC) sentences ended with a semantically anomalous word which was totally unexpected given the sentential context, e.g., _De klimmers bereikten eindelijk de top van de tulp_ (_The climbers finally reached the top of the tulip_). The third type of sentences ended with a semantically anomalous word that had the same initial phonemes (and lexical stress) as the high-cloze words from the congruent condition: initially congruent (IC). There were 87 trials per condition for each of the three conditions, and a set of 87 filler sentences were added. From the EEG literature it is known that a stronger negative potential is produced following incongruent compared to congruent sentence endings about 300-500 ms after the word onset. This response is termed the N400 effect¹ ². For more information about the materials you could take a look at the published EEG experiment using the same sentence materials³. + +In the study applied here, the subjects were seated in a relaxed position under the MEG helmet. Their task was to attentively listen to spoken sentences. They were informed that some of the sentences would be semantically anomalous. Acoustic transducers were used to deliver the auditory stimuli. After a 300-ms warning tone, followed by a 1200 ms pause, a sentence was presented. Every next trial began 4100 ms after the offset of the previous sentence. To reduce eye blinks and movements in the time interval in which the sentence was presented, subjects were instructed to fixate on an asterisk presented visually 1000 ms prior to the beginning of the sentence. The asterisk remained on the screen until 1600 ms after the onset of the spoken sentence. Subjects were encouraged to blink when the asterisk was not displayed on the screen. + +MEG signals were recorded with a 151-channel CTF system. 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 100 Hz, digitized at 300 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. + +The MEG data are stored as epochs or trials of fixed length around each stimulus trigger. + +## Background on time-frequency analysis + +Oscillatory components contained in the ongoing EEG or MEG signal often show power changes relative to experimental events. These signals are not necessarily phase-locked to the event and will not be represented in event-related fields and potentials ([Tallon-Baudry and Bertrand (1999)](https://doi.org/10.1016/S1364-6613(99)01299-1)). The goal of this section is to compute and visualize event-related changes by calculating time-frequency representations (TFRs) of power. This will be done using analysis based on Fourier analysis and wavelets. The Fourier analysis will include the application of multitapers ([Mitra and Pesaran (1999)](https://doi.org/10.1016/S0006-3495(99)77236-X), [Percival and Walden (1993)](http://lccn.loc.gov/92045862)) which allow a better control of time and frequency smoothing. + +Calculating time-frequency representations of power is done using a sliding time window. This can be done according to two principles: either the time window has a fixed length independent of frequency, or the time window decreases in length with increased frequency. For each time window the power is calculated. Prior to calculating the power one or more tapers are multiplied with the data. The aim of the tapers is to reduce spectral leakage and control the frequency smoothing. + +{% include image src="/assets/img/workshop/neuroimaging2-2425/sped4/figure1.png" width="600" %} + +_Figure: Time and frequency smoothing. (a) For a fixed length time window the time and frequency smoothing remains fixed. (b) For time windows that decrease with frequency, the temporal smoothing decreases and the frequency smoothing increases._ + +If you want to know more about tapers/ window functions you can have a look at this +[Wikipedia site](https://en.wikipedia.org/wiki/Window_function). Note that Hann window is another name for Hanning window used in this tutorial. There is also a Wikipedia site about multitapers, to take a look at it click [here](https://en.wikipedia.org/wiki/Multitaper). + +## Procedure + +To calculate the time-frequency analysis for the example dataset we will perform the following steps: + +- Read the data into MATLAB using **[ft_definetrial](/reference/ft_definetrial)** and **[ft_preprocessing](/reference/ft_preprocessing)** +- Seperate the trials from each condition using **[ft_selectdata](/reference/utilities/ft_selectdata)** +- Compute the power values for each frequency bin and each time bin using the function **[ft_freqanalysis](/reference/ft_freqanalysis)** +- Visualize the results. This can be done by creating time-frequency plots for one (**[ft_singleplotTFR](/reference/ft_singleplotTFR)**) or several channels (**[ft_multiplotTFR](/reference/ft_multiplotTFR)**), or by creating a topographic plot for a specified time- and frequency interval (**[ft_topoplotTFR](/reference/ft_topoplotTFR)**). + +{% include image src="/assets/img/workshop/neuroimaging2-2425/sped4/figure2.png" width="200" %} + +_Figure: Schematic overview of the steps in time-frequency analysis_ + +In this tutorial, procedures of 4 types of time-frequency analysis will be shown. You can see each of them under the titles: Time-frequency analysis I, II ... and so on. + +## Preprocessing + +The first step is to read the data using the function **[ft_preprocessing](/reference/ft_preprocessing)**. It is recommended to read larger time intervals than the time period of interest. In this example, the time of interest is from -0.5 s to 1.5 s (t = 0 s defines the time of stimulus); however, the script reads the data from -1.0 s to 2.0 s. + +{% include markup/skyblue %} +**Exercise 1**: Why is it recommended to read in larger time intervals than the time window of interest when we want to do time-frequency analysis? +{% include markup/end%} + +### Reading in the data + +First, [download the data from here](https://download.fieldtriptoolbox.org/tutorial/Subject01.zip) and unzip it somewhere. Make sure FieldTrip is added to your MATLAB path: + + addpath + ft_defaults + +Then, execute the following code, which will determine the time indices of the trials of interest (note: adapt the path to wherever you unzipped Subject01.ds): + + cfg = []; + cfg.dataset = '/Subject01.ds'; + cfg.trialfun = 'ft_trialfun_general'; % this is the default + cfg.trialdef.eventtype = 'backpanel trigger'; % name of the trigger channel in the dataset + cfg.trialdef.eventvalue = [3 5 9]; % the values of the stimulus trigger for the three conditions + % 3 = fully incongruent (FIC), 5 = initially congruent (IC), 9 = fully congruent (FC) + cfg.trialdef.prestim = 1; % in seconds + cfg.trialdef.poststim = 2; % in seconds + + cfg = ft_definetrial(cfg); + +{% include markup/skyblue %} +**Exercise 2**: What do the configuration values `cfg.trialdef.prestim` and `cfg.trialdef.poststim` denote? Have a look at the documentation by doing `edit ft_definetrial` (or browse this wiki) if needed. +{% include markup/end%} + +The resulting `cfg` structure will have a `cfg.trl` matrix which contains the begin and end samples of all trials of interest, as well as the offset (3rd column) that determines which sample corresponds to time points t = 0. + +### Cleaning + +Some trials have previously been identified as artifactual (due to for example eye blinks or MEG SQUID jumps). Also, two MEG channels were malfunctioning. Both these trials and channels need to be removed. Furthermore, while reading in the data, we remove the overall per-trial and per-channel mean to facilitate downstream time-frequency analysis. The following code achieves all this and reads in the data: + + % remove the trials that have artifacts from the trl + cfg.trl([2, 5, 6, 8, 9, 10, 12, 39, 43, 46, 49, 52, 58, 84, 102, 107, 114, 115, 116, 119, 121, 123, 126, 127, 128, 133, 137, 143, 144, 147, 149, 158, 181, 229, 230, 233, 241, 243, 245, 250, 254, 260],:) = []; + + % 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 + + data_all = ft_preprocessing(cfg); + +For subsequent analysis in this hands-on, we extract the trials of the fully incongruent condition: + + cfg = []; + cfg.trials = data_all.trialinfo == 3; + dataFIC = ft_redefinetrial(cfg, data_all); + +Subsequently you can save the data to disk, to easily continue later on, without having to read in all data again: + + save dataFIC dataFIC + +## TFR I: Hanning taper, fixed window length + +Here, we will describe how to calculate time frequency representations using the commonly used [Hanning, or Hann, tapers](https://en.wikipedia.org/wiki/Hann_function). When choosing for a fixed window length procedure, the frequency resolution is defined according to the length of the time window (delta T). The frequency resolution (delta f in figure 1) = 1/length of time window in sec (delta T in figure 1). Thus a 500 ms time window results in a 2 Hz frequency resolution (1/0.5 sec = 2 Hz) meaning that power can be calculated for 2 Hz, 4 Hz, 6 Hz etc. An integer number of cycles should fit in the time window. + +In the following example, a time window with length 500 ms is applied: + + cfg = []; + cfg.output = 'pow'; + cfg.channel = 'MEG'; + cfg.method = 'mtmconvol'; + cfg.taper = 'hanning'; + cfg.foi = 2:2:30; % analysis 2 to 30 Hz in steps of 2 Hz + cfg.t_ftimwin = ones(length(cfg.foi),1).*0.5; % length of time window = 0.5 sec + cfg.toi = -0.5:0.05:1.5; % time window "slides" from -0.5 to 1.5 sec in steps of 0.05 sec (50 ms) + TFRhann = ft_freqanalysis(cfg, dataFIC); + +The field `cfg.method = 'mtmconvol';` instructs `ft_freqanalysis` to execute a convolution-type analysis. (The "mtm" stands for "multitaper method", which will be explained in more detail below, but is mainly there for historical reasons, since the function also supports convolution that is *not* based on multitapers.) Regardless of the method used for calculating the TFR, the output format is identical. It is a structure with the following fields: + + TFRhann = + label: {149x1 cell} % Channel names + dimord: 'chan_freq_time' % Dimensions contained in powspctrm, channels X frequencies X time + freq: [2 4 6 8 10 12 14 16 18 20 22 24 26 28 30] % Array of frequencies of interest (the elements of freq may be different from your cfg.foi input depending on your trial length) + time: [1x41 double] % Array of time points considered + powspctrm: [149x15x41 double] % 3-D matrix containing the power values + elec: [1x1 struct] % Electrode positions etc + grad: [1x1 struct] % Gradiometer positions etc + cfg: [1x1 struct] % Settings used in computing this frequency decomposition + +The field `TFRhann.powspctrm` contains the temporal evolution of the raw power values for each specified frequency. + +## Visualization with standard MATLAB code + +To facilitate understanding of the output of `ft_freqanalysis`, it is instructive to first plot it all in the same way we did before using `imagesc`. To do this, we need to first average over channels to obtain a 2D matrix. + + pow_allchan = squeeze(mean(TFRhann.powspctrm, 1)); + figure; + imagesc(TFRhann.time, TFRhann.freq, pow_allchan); + axis xy; + colorbar(); + xlabel('time (s)'); + ylabel('frequency (Hz)'); + +{% include markup/skyblue %} +**Exercise 3**: What do you notice most clearly in these raw power values? How does power vary across frequencies? And across time? +{% include markup/end%} + +## Visualization with FieldTrip code + +As discussed in the lectures, biological signals are often dominated by a strong so-called "1/f" component. In order to visualize and interpret task-related *changes* in oscillatory activity with respect to a baseline window, it is therefore recommended to perform baseline normalization. + +In general there are two possibilities for normalizing: + +- Subtracting, for each frequency, the average power in a baseline interval from all other power values. This gives, for each frequency, the absolute change in power with respect to the baseline interval. +- Expressing, for each frequency, the raw power values as the relative increase or decrease with respect to the power in the baseline interval. This means active period/baseline. If we furthermore log-transform (and scale) these ratios, we end up with the commonly used decibel (db). + +There are three ways of graphically representing the data: 1) time-frequency plots of all channels, in a quasi-topographical layout, 2) time-frequency plot of an individual channel (or average of several channels), 3) topographical 2-D map of the power changes in a specified time-frequency interval. + +To plot the TFRs from all the sensors use the function **[ft_multiplotTFR](/reference/ft_multiplotTFR)**. Settings can be adjusted in the cfg structure. For example: + + cfg = []; + cfg.baseline = [-0.5 -0.1]; + cfg.baselinetype = 'db'; % use db baseline correction + cfg.zlim = 'maxabs'; % color scale symmetric around zero + cfg.showlabels = 'yes'; + cfg.layout = 'CTF151_helmet.mat'; + figure + ft_multiplotTFR(cfg, TFRhann); + +Note that using the options `cfg.baseline` and `cfg.baselinetype` results in baseline correction of the data, implicitly during the plotting call. Baseline correction can also be applied explicitly by calling **[ft_freqbaseline](/reference/ft_freqbaseline)**. Moreover, you can combine the various visualization options/functions interactively to explore your data. + +An interesting effect seems to be present in the TFR of sensor MRC15. To make a plot of a single channel use the function **[ft_singleplotTFR](/reference/ft_singleplotTFR)**. + + cfg = []; + cfg.baseline = [-0.5 -0.1]; + cfg.baselinetype = 'db'; + cfg.zlim = 'maxabs'; + cfg.channel = 'MRC15'; + cfg.layout = 'CTF151_helmet.mat'; + figure + ft_singleplotTFR(cfg, TFRhann); + +From the previous figure you can see that there is an increase in power around 15-20 Hz in the time interval 0.9 to 1.3 s after stimulus onset. To show the topography of the beta increase use the function **[ft_topoplotTFR](/reference/ft_topoplotTFR)**. + + cfg = []; + cfg.baseline = [-0.5 -0.1]; + cfg.baselinetype = 'db'; + cfg.xlim = [0.9 1.3]; + cfg.zlim = 'maxabs'; + cfg.ylim = [15 20]; + cfg.marker = 'on'; + cfg.layout = 'CTF151_helmet.mat'; + cfg.colorbar = 'yes'; + figure + ft_topoplotTFR(cfg, TFRhann); + +{% include markup/skyblue %} +**Exercise 4**: By default, FieldTrip plotting functions support an interactive mode. (Note: This interactive mode does not work in combination with Matlab "Live scripts", if you happen to be using these.) The interactive mode allows you to drag a box around sensors of interest, click that box, and you'll get an average TFR for those sensors only. In the resulting TFR plot, drag a box around a time/frequency window of interest, and you'll see a topographical plot. From that topoplot, you can again select sensors and go to an averaged TFR, etc. Optionally see also the [plotting tutorial](/tutorial/plotting) for more details. Play around with interactive mode after `ft_multiplotTFR` and reflect briefly on what you see. +{% include markup/end%} + +{% include markup/skyblue %} +**Exercise 5**: Plot the TFR of sensor MLC24. How do you account for the increased power at ~300 ms post-stimulus (hint: compare to what you might expect in an event-related field)? +{% include markup/end%} + +## TFR II: Hanning taper, frequency dependent window length + +It is also possible to calculate the TFRs with respect to a time window that varies with frequency. Typically the time window gets shorter with an increase in frequency. The main advantage of this approach is that the temporal smoothing decreases with higher frequencies, leading to increased sensitivity to short-lived effects. However, an increased temporal resolution is at the expense of frequency resolution. + +{% include markup/skyblue %} +**Exercise 6**: Why is it the case that increased temporal resolution is at the expense of frequency resolution? +{% include markup/end%} + +We will here show how to perform a frequency-dependent time-window analysis, using a sliding window Hanning taper based approach. The approach is very similar to wavelet analysis. A wavelet analysis performed with a Morlet wavelet mainly differs by applying a Gaussian shaped taper (see [Time-frequency analysis IV](#time-frequency-analysis-iv)). + +The analysis is best done by first selecting the numbers of cycles per time window which will be the same for all frequencies. For instance if the number of cycles per window is 7, the time window is 1000 ms for 7 Hz (1/7 x 7 cycles); 700 ms for 10 Hz (1/10 x 7 cycles) and 350 ms for 20 Hz (1/20 x 7 cycles). The frequency can be chosen arbitrarily - however; too fine a frequency resolution is just going to increase the redundancy rather than providing new information. + +Below is the configuration for a 7-cycle time window. The calculation is only done for one sensor (MRC15) but it can of course be extended to all sensors. + + cfg = []; + cfg.output = 'pow'; + cfg.channel = 'MRC15'; + cfg.method = 'mtmconvol'; + cfg.taper = 'hanning'; + cfg.foi = 2:1:30; + cfg.t_ftimwin = 7./cfg.foi; % 7 cycles per time window + cfg.toi = -0.5:0.05:1.5; + TFRhann7 = ft_freqanalysis(cfg, dataFIC); + +To plot the result use **[ft_singleplotTFR](/reference/ft_singleplotTFR)**: + + cfg = []; + cfg.baseline = [-0.5 -0.1]; + cfg.baselinetype = 'db'; + cfg.zlim = 'maxabs'; + cfg.channel = 'MRC15'; + cfg.interactive = 'no'; + cfg.layout = 'CTF151_helmet.mat'; + figure + ft_singleplotTFR(cfg, TFRhann7); + +{% include markup/skyblue %} +**Exercise 7**: Adjust the length of the time-window and thereby degree of smoothing. Use **[ft_singleplotTFR](/reference/ft_singleplotTFR)** to show the results. Discuss the consequences of changing these setting. +{% include markup/end%} + +## TFR III: Morlet wavelets + +As discussed in detail in the lectures, a common way to calculate TFRs is convolution with Morlet wavelets. The approach is equivalent to calculating TFRs with sliding time windows that depend on frequency using a taper with a Gaussian shape. + +{% include markup/skyblue %} +**Exercise 8**: Why are the two approaches equivalent? (Approach 1: slide a time window over your data, multiply the data in each window with a Gaussian, then FFT; approach 2: construct a wavelet by multiplying a complex sinusoid with a Gaussian window, and convolve that wavelet with your data.) +{% include markup/end%} + +The commands below illustrate how to do Morlet-wavelet based analysis in FieldTrip. One crucial parameter to set is `cfg.width`. It determines the width of the wavelets in number of cycles. Making the value smaller will increase the temporal resolution at the expense of frequency resolution and vice versa. + +Calculate TFRs using Morlet wavelet convolution: + + cfg = []; + cfg.channel = 'MEG'; + cfg.method = 'wavelet'; + cfg.width = 7; + cfg.output = 'pow'; + cfg.foi = 1:2:30; + cfg.toi = -0.5:0.05:1.5; + TFRwave = ft_freqanalysis(cfg, dataFIC); + +Plot the result: + + cfg = []; + cfg.baseline = [-0.5 -0.1]; + cfg.baselinetype = 'db'; + cfg.zlim = 'maxabs'; + cfg.showlabels = 'yes'; + cfg.layout = 'CTF151_helmet.mat'; + cfg.colorbar = 'yes'; + figure + ft_multiplotTFR(cfg, TFRwave) + +{% include markup/skyblue %} +**Exercise 9A**: Adjust `cfg.width` and see how the TFRs change. **Exercise 9B**: Make some plots using `absolute` baseline correction instead of `db` and see how the TFRs change. I'd recommend switching to `ft_singleplotTFR` for one or a few channels of interest for these exercises, rather than doing the full `ft_multiplotTFR` each time. +{% include markup/end%} + +## TFR IV: Multitapers + +Multitapers (literally: multiple tapers per time window of interest) are sometimes used in order to achieve better control over the frequency smoothing. More tapers for a given time window will result in stronger smoothing. For frequencies above 30 Hz, smoothing has been shown to be advantageous, increasing sensitivity thanks to reduced variance in the estimates despite reduced effective spectral resolution. Oscillatory gamma activity (30-100 Hz) is quite broad band and thus analysis of this signal component benefits from multitapering, which trades spectral resolution against increased sensitivity. For signals lower than 30 Hz it is recommend to use only a single taper, e.g., a Hanning taper as shown above (beware that in the example below multitapers are used to analyze low frequencies because there are no effects in the gamma band in this dataset). + +Time-frequency analysis based on multitapers is also performed by **[ft_freqanalysis](/reference/ft_freqanalysis)**. The function uses a sliding time window for which the power is calculated for a given frequency. Prior to calculating the power by discrete Fourier transforms the data are 'tapered'. Several orthogonal tapers might be used for each time window. The power is calculated for each tapered data segment and then combined. In the example below we apply a time window which gradually becomes shorter for higher frequencies (similar to wavelet techniques). Note that this is not necessary, but up to the researcher to decide. The arguments for the chosen parameters are as follows: + +- `cfg.foi`, the frequencies of interest, here from 1 Hz to 30 Hz in steps of 2 Hz. The step size could be decreased (smaller steps) at the expense of computation time and redundancy, but leading to smoother plots. +- `cfg.toi`, the time-interval of interest. This vector determines the center times for the time windows for which the power values should be calculated. The setting `cfg.toi = -0.5:0.05:1.5` results in power values from -0.5 to 1.5 s in steps of 50 ms. A finer time resolution will give redundant information and longer computation times, but a smoother graphical output. +- `cfg.t_ftimwin` is the length of the sliding time-window in seconds. We have chosen `cfg.t_ftimwin = 5./cfg.foi`, i.e. 5 cycles per time-window. When choosing this parameter it is important that a full number of cycles fit within the time-window for a given frequency. +- `cfg.tapsmofrq` determines the width of frequency smoothing in Hz. We have chosen `cfg.tapsmofrq = cfg.foi*0.4`, i.e. the smoothing will increase with frequency. Specifying larger values will result in more frequency smoothing. For less smoothing you can specify smaller values. + +These settings result in the following characteristics as a function of the frequencies of interest: + +{% include image src="/assets/img/workshop/neuroimaging2-2425/sped4/figure7.png" width="400" %} + +_Figure: a) The characteristics of the TFRs settings using multitapers in terms of time and frequency resolution of the settings applied in the example. b) Examples of the time-frequency tiles resulting from the settings._ + + cfg = []; + cfg.output = 'pow'; + cfg.channel = 'MEG'; + cfg.method = 'mtmconvol'; + cfg.foi = 1:2:30; + cfg.t_ftimwin = 5./cfg.foi; + cfg.tapsmofrq = 0.4 *cfg.foi; + cfg.toi = -0.5:0.05:1.5; + TFRmult = ft_freqanalysis(cfg, dataFIC); + +Note that if you get an error here that `dpss` is not working because you require the Signal Processing Toolbox, you can either (a) install the Signal Processing Toolbox if you have a license for it, or manually add `/external/signal/dpss_hack` to your Matlab path. + +Plot the result (again in the command window for interactive plotting): + + cfg = []; + cfg.baseline = [-0.5 -0.1]; + cfg.baselinetype = 'db'; + cfg.zlim = 'maxabs'; + cfg.showlabels = 'yes'; + cfg.layout = 'CTF151_helmet.mat'; + cfg.colorbar = 'yes'; + figure + ft_multiplotTFR(cfg, TFRmult) + +{% include markup/skyblue %} +**Exercise 10**: Explore the TFRs that result from multitapering. Reflect on how the alpha/beta-band (~10-20 Hz) response around 1-1.5s post-stimulus appears now, in comparison with the earlier approaches. Why might it be beneficial to use multitapering? +{% include markup/end%} + +### Multitapering as a hack around the time-frequency uncertainty principle? + +A final more detailed note on what multitapering actually does. While the fundamental trade-off between time- and frequency-resolution cannot be broken, multitapering offers a sort of "hack" to artificially reduce one *without* increasing the other. Sometimes we want to smooth over different frequencies, which is another way of saying that sometimes we *want* to reduce our frequency resolution. For example, we might know that, from a cognitive/physiological perspective, the exact same phenomenon is reflected in slightly different frequency bands across participants. By now you may have the (correct) intuition that in order to do this, we should reduce our time window or wavelet width (since that increases time resolution and thereby decreases frequency resolution a.k.a. increases frequency smoothing). However, we may also know, again from a cognitive/physiological perspective, that the exact same phenomenon is not always present at the exact same time points across participants, or even across trials! So here, actually, increasing our temporal resolution (a.k.a. decreasing time smoothing) is not what we want, because also that reduces our subsequent statistical sensitivity. Multitapering offers a way to reduce our frequency resolution (increase smoothing) while keeping the time window the same. Note that this is not a magical way around the "uncertainty principle", as it can only *increase* smoothing beyond that which is inherent in the time window, it cannot decrease it (i.e., increase frequency resolution) beyond that which is dictated by the time window length. + + + + +## References + +1. Kutas M, Hillyard SA. (1980) Reading senseless sentences: brain potentials reflect semantic incongruity. Science. 207(4427):203-5 +2. Kutas M, Federmeier KD. (2000) Electrophysiology reveals semantic memory use in language comprehension. Trends Cogn Sci. 4(12):463-470 +3. van den Brink D, Brown CM, & Hagoort P. (2001). Electrophysiological evidence for early contextual influences during spoken-word recognition: N200 versus N400 effects. J Cogn Neurosci. 13(7):967-985 +4. Wang L, Jensen O, van den Brink D, Weder N, Schoffelen JM, Magyari L, Hagoort P, Bastiaansen M. (2012) Beta oscillations relate to the N400m during language comprehension. Hum Brain Mapp. 2012 Dec;33(12):2898-912. diff --git a/workshop/neuroimaging2-2526/sr1.md b/workshop/neuroimaging2-2526/sr1.md new file mode 100644 index 000000000..707a86901 --- /dev/null +++ b/workshop/neuroimaging2-2526/sr1.md @@ -0,0 +1,339 @@ +--- +title: SR1 - Source reconstruction +tags: [neuroimaging2-2425] +--- + +## 1 The biophysical forward model and linear mixing + +We assume that you have MATLAB installed and that you understand the basics +of it. This includes the following: starting a MATLAB session, using the command-line, +and using the MATLAB-editor. + +For the exercises we are going to use MATLAB-code that has been written specifically +for this course, as well the MATLAB-based toolbox FieldTrip. The latter has +been initiated at the Donders Centre for Cognitive Neuroimaging. The toolbox +is under continuous development and contains many state-of-the-art algorithms +for the analysis of electrophysiological data. Unlike other popular toolboxes, +FieldTrip does not come with a Graphical User Interface (GUI), and thus simply +clicking buttons is not sufficient to achieve what comes to your mind. As a +consequence, the novice user may experience some difficulties in getting the +most out of the toolbox. Fortunately, there is an excellent website that provides +a wealth of information. For the exercises in this course, no detailed knowledge +about FieldTrip is required. However, if you have some time to spare and feel +adventurous, there is a large quantity of [tutorials](/tutorial) that can be followed. + +### 1.1 Getting started: setting up the MATLAB environment + +As a first step, you need to install the required toolboxes on your computer. +These can be obtained from Brightspace or GitHub. We advise you to create a +separate folder that will contain the material for this course. For example, +on Windows, you could create a folder in 'My Documents' called 'neuroimaging2', +with a subfolder called 'matlab'. You have to download the two zip-files `fieldtrip.zip` +and `ni2.zip` from Brightspace to this folder and unzip them in this location. + +Next, MATLAB needs to know where to find the code we are using in the exercises. +To this end, the folders you have just created need to be added to the MATLAB-path. +Although this can be done through the 'Set Path' option (see figure 1), we will +use a `startup.m` file. This is a script that automatically deals with the path-settings +that are relevant for the exercises. Each time you start a MATLAB-session, you +need to change the Current Folder to `ni2`, and type `ni2_startup` on the command +line. This will automatically ensure that the path-settings are correct. + + ni2_startup + +### 1.2 Getting started: tips on how to approach the exercises + +The best way to approach the exercises is to create a MATLAB-script that can +be used to copy-and-paste the MATLAB-code in the instructions, and allows you +to play around with various parameter settings etc. You can open a new script +in the MATLAB-editor by typing edit on the command line, or by going to File->New->Script. +Next, you can save each script as a separate file in your course folder. + +## 2 Spatial mixing is linear mixing is matrix multiplication + +Electrophysiological signals typically represent a mixture of the underlying +neuronal sources. The next exercises intend to give some fundamental insights +with respect to the general mixing phenomenon. After doing these exercises: + +- You will understand that instantaneous mixing is the same as making a linear combination of the underlying 'source' processes. +- You will know that making a weighted combination is equivalent to doing a matrix multiplication. +- Argue that temporal characteristics of the sources and the spatial mixing interact in a non-trivial way. + +We start with the generation of some 'source signals'. To this end we can +use the function `ni2_activation`. If you type `help ni2_activation` on the +MATLAB-command line you can read how to use the function. Let's start with + + [sourceactivity, time] = ni2_activation; + whos sourceactivity time + +Without specifying any input arguments, the function will set all parameters +to a default value. Type `help ni2_activation` to see the default parameter +values. + +The output variables `sourceactivity` and `time` contain for each sampled +time point the amplitude of the source signal, and the time in seconds. You +can visualize the activity by using MATLAB's built-in function plot. Type + + figure; plot(time, sourceactivity); xlabel('time (s)'); ylabel('activation (a.u.)') + +Imagine that a stimulus was presented at time t=0, and that around 300 or +400 ms after that stimulus some cortical activation was triggered. + +**_Verify whether the default parameter settings in ni2_activation coincide +with what is described in the function's._** + +Now we will create a second time course that peaks at 0.4 seconds and which +has an oscillation frequency of 5 Hz: + + [sourceactivity2, time] = ni2_activation('latency', 0.4, 'frequency', 5); + + whos sourceactivity2 time + +Verify the result by using the plot function. You can visualize the two time +courses in the same figure by doing: + + figure; hold on; xlabel('time (s)'); ylabel('activation (a.u.)') + plot(time, sourceactivity, 'b'); + plot(time, sourceactivity2, 'r'); + +The `hold on` command is necessary in order to show more than one line in +the figure. + +**_Evaluate the effect of the different parameter settings is on the shape +of the time course, by specifying different values for `frequency`, `ncycle`, +`latency` and `phase`._** + +We can apply a very simple instantaneous mixing to the two sources' time courses +by just summing the two. This procedure is equivalent to taking a weighted combination, +where the weight for each source is 1: + + datamix = sourceactivity + sourceactivity2; + size(datamix) + figure; plot(time, datamix); + +**_Plot the three time courses in a single figure and verify that the peaks +and troughs in the mixed signal do not necessarily end up at the same latencies +as in the constituent source signals._** + +**_Investigate how the morphology of the mixed signal changes, when the latency +of the second component (`sourceactivity2`) is changed. Also, see what the effect +is of changing the phase of the second component. Hint: create a new `sourceactivity2` +with different input parameters to `ni2_activation` and create a new `datamix`. +Look in the help documentation of `ni2_activation` to try and understand how to +adjust the phase of the signal._** + +Next, we will verify that a (weighted) combination of two (or more) activation +time courses can be easily achieved by means of a matrix multiplication. To +this end we first concatenate the original source activations into a single +matrix (each of the original time courses will be a row in the resulting matrix): + +make sure that the `sourceactivity2` is what it should be, not the one you have been playing with + + [sourceactivity2, time] = ni2_activation('latency', 0.4, 'frequency', 5); + +concatenate the two vectors in a matrix + + sourcecombined = [sourceactivity; sourceactivity2]; + size(sourcecombined) + +Note that it is now very easy to plot both time courses with a single plotting +command: + + figure + plot(time, sourcecombined); + +We can now make a mixture of the two 'sources' with a simple matrix multiplication. +Execute the following lines of code and verify that the outcomes are the same. + + datamix1 = sourceactivity + sourceactivity2; + datamix2 = sum(sourcecombined, 1); + datamix3 = [1 1] * sourcecombined; + + whos datamix1 datamix2 datamix3 + +or, more generally + + mix = [1 1]; + datamix3 = mix * sourcecombined; + +Note that the outcome of a matrix multiplication A*B, where the input matrices +are of dimensionality n x m and m x p yields a matrix of dimensionality n x +p. If you don't remember exactly the mechanics involved in a matrix multiplication, +you should look this up in the course material for the CNS Advanced Mathematics +course, or you could check the Wikipedia entry on matrix multiplication. In +our example p corresponds with the number of time points, m with the number +of 'sources', and n with the number of 'channels'. Our (very simple) 1 x 2 mixing +matrix [1 1] thus corresponds to the mapping of 2 sources onto a single channel, +where each of the sources is weighted with a factor of 1. We can increase the +number of channels by increasing the number of rows in the matrix mix: + + mix = [0 1; 0.1 0.9; 0.25 0.75; 0.5 0.5; 0.75 0.25; 0.9 0.1; 1 0]; + datamix = mix * sourcecombined; + whos mix + whos sourcecombined + whos datamix + +To better evaluate the individual time courses of the mixed sources, but still +visualizing them in the same figure, you can type: + + figure; plot(time, datamix + repmat((1:7)', [1 1000])); + +The part _repmat((1:7)', [1 1000]))_ lifts the 7 signals to different heights +in the figure. To better see what the repmat function does, type + + figure; plot(time, repmat((1:7)', [1 1000])); + +Note that the sum of the values in the rows of the mixing matrix do not have +to sum to a value of one. Note also that the values in the mixing matrix can +be negative. + +**_What is the effect of a negative mixing weight on the mixed result?_** + +## 3 Spatial mixing of neural sources: the leadfield + +This section deals with the spatial topographies of activity that are observed +with EEG when neural sources are active. Modeling these spatial topographies +when the underlying sources are known is called forward modeling. After doing +these exercises + +- You will understand the concept of a leadfield. +- You understand the effect of source parameters on the observed topographies. + +To build a forward model one needs a model of the sources, a volume conduction +model, and information about the sensor array (we will come back to this somewhat +later in this course). Sources are typically modeled as equivalent current dipoles +(ECDs). At the spatial scale of EEG/MEG an ECD can be represented as a point +source with 6 parameters: 3 location parameters (x, y, z coordinates in a (Cartesian) +coordinate system), and 3 parameter dipole moment (combining amplitude and orientation +information) along the x, y, z-axes of the coordinate system (see figure). + +As a volume conduction model we will initially use a simple model, consisting +of three concentric spheres. These spheres represent the inner and outer surfaces +of the skull, and the scalp surface. In an EEG volume conduction model, these +surfaces represent the locations at which the conductivity of the tissue changes +rapidly, leading to distortion and smearing of the volume currents. We will +use an EEG sensor array that has 91 electrodes placed on the spherical scalp +surface. As such this sensor array does not correspond to what you typically +encounter in real life EEG-setups, but this is not a problem for bringing across +the main message. + +To create the headmodel, type: + + headmodel = ni2_headmodel('type', 'spherical', 'nshell', 3) + +To create the sensor array, type: + + sensors = ni2_sensors('type', 'eeg') + +These can be visualized using a few functions from the FieldTrip-toolbox: + + figure; hold on; + ft_plot_headmodel(headmodel, 'edgecolor', 'none', 'facecolor', 'skin'); alpha 0.5 + ft_plot_sens(sensors, 'elecshape', 'disc', 'elecsize', 1, 'label', 'on'); + ft_plot_axes(sensors); + view([1 1 1]) + axis on; grid on + xlabel('x'); ylabel('y'); zlabel('z'); + +Now that we have a sensor array and a volume conductor model of the head, +we can build a model for any source specified within the volume conductor. For +example: + + dippar1 = [0 0 6 1 0 0]; % dipole position at [0 0 6], dipole moment is [1 0 0] + leadfield1 = ni2_leadfield(sensors, headmodel, dippar1); + ni2_topoplot(sensors, leadfield1); colorbar + +The variable dippar1 consists of 6 elements, where the first triplet of numbers +represents the x, y, z-position of the dipole, and the second triplet of numbers +represents the dipole moment. + +**_Take a few moments to try and understand the topography._** + +A basic feature of a leadfield is that a change in the amplitude of the dipole +does not change the distribution of the peaks in the topography. This can be +explored as follows: + + dippar2 = [0 0 6 2 0 0]; + leadfield2 = ni2_leadfield(sensors, headmodel, dippar2); + ni2_topoplot(sensors, leadfield2); colorbar + +**_Verify that the values in leadfield2 differ from the values in leadfield1 +by a factor of 2._** + +The important consequence of this is that the amplitude can be 'uncoupled' +from the orientation (and location) information. Namely: + + ni2_topoplot(sensors, leadfield1*2); colorbar + +will yield the same result as: + + ni2_topoplot(sensors, leadfield2); colorbar + +For this reason, it is custom (and convenient) that leadfields are defined +as the spatial topography of the electric potential (EEG) or magnetic field +(MEG) due to a unit-amplitude source at a given location and orientation. + +In the example above the amplitude information (the value '2') just represented +the activity of a source in a single time point. The amplitude term can be very +easily extended to contain multiple time points, i.e., to become a vector. Mathematically, +this can still be expressed as a matrix multiplication: + + leadfield * amplitude + +We will return to this very important property in the next section. + +Another basic feature is that if two (or more sources) are simultaneously +active, the topographies caused by the individual sources mix linearly to yield +a composite topography: + + dippar3 = [0 0 5.99 1 0 0; 0 0 6.01 1 0 0]; % two dipoles at almost exactly the same location and with the same moment + leadfield3 = ni2_leadfield(sensors, headmodel, dippar3); + topo = leadfield3*[1; 1]; + ni2_topoplot(sensors, topo); colorbar + +The variable _dippar3_ consists of 2 rows of 6 parameters, each representing +a dipole. As a result, the variable leadfield3 consists of two columns, each +representing the spatial topography that is caused by one dipole. + +Note that the mixing column vector ([1; 1]) features on the right side of +the leadfield matrix, as opposed to on the left side (where in addition the +mixing vector was a row vector) when we were mixing activation time courses +in the previous section. Note also, that in the above example two unit-amplitude +sources are located very close to one another. + +**_Verify that the topography due to 2 very closely spaced sources (with the +same orientation) is hardly distinguishable from a situation where only one +source (with twice the single sources’ amplitude) is present at the average +location of the two sources._** + +This is a very important feature of electromagnetic forward models and in +essence the cause for the fact that the inverse problem is ill-posed: there +is no unique solution. We will revisit this issue in more detail at a later +stage. + +So far, we have considered the leadfield of a source located at [0 0 6], with +an orientation of [1 0 0]. + +**_Explore the effect of changing the position and orientation parameters. +Investigate the following situations:_** + +1. **_Keep the orientation fixed and change the location of the source, e.g., +by moving it closer to the center of the 'head’._** +2. **_Change the y and z position of the source, keeping x the same._** +3. **_Keep the location at a fixed position and change the orientation of the +source. Start by investigating the orientations [1 0 0], [0 1 0], and [0 0 1]._** +4. **_Also explore source orientations that have more than 1 non-zero value, +e.g. [0 1 1]._** +5. **_Compare the effect of a change in the sign of the orientation, e.g. [0 +1 1] versus [0 -1 -1]._** + +Up to this point we have explicitly added the orientation parameters when +computing a leadfield. Typically, however, a leadfield is defined as a 3-column +matrix, where, for a source at a given location, each of the columns represent +a unit source with the orientation along the axes of the coordinate system. +In other words, the first column represents a source with orientation [1 0 0], +the second column represents a source with orientation [0 1 0], etc. Based on +this 3-column representation, any orientation of the source can be represented +by means of (could it be any different?) linear mixing. + +**_Verify this._** diff --git a/workshop/neuroimaging2-2526/sr2.md b/workshop/neuroimaging2-2526/sr2.md new file mode 100644 index 000000000..8a0222c17 --- /dev/null +++ b/workshop/neuroimaging2-2526/sr2.md @@ -0,0 +1,333 @@ +--- +title: SR2 - Source reconstruction +tags: [neuroimaging2-2425] +--- + +## 4 Spatiotemporal mixing of neural sources: the basic observation model + +This document continue with the MATLAB exercises that form part of the course +“Neuroimaging II” relating to forward modeling. + +The best way to approach the exercises is to create a MATLAB-script that can +be used to copy-and-paste the MATLAB-code in the instructions, and allows you +to play around with various parameter settings etc. You can open a new script +in the MATLAB-editor by typing edit on the command line, or by going to File->New->Script. +Next, you can save each script as a separate file in your course folder. + +Set up the MATLAB path + + ni2_startup + ft_version + +_This section deals with the basic observation model, and combines what we +have learnt in the previous section. The basic observation model relates to +the very generic model of how the observed data are generated by the constituent +sources. It can be briefly stated as `y = L*s + n`. Here, `y` is the sensor +data, `L` is the known leadfield, `s` is the source amplitude and `n` is sensor +noise._ + +_After doing these exercises:_ + +- _You understand this very basic equation._ +- _In particular, you will understand that the sources’ time courses can be +'uncoupled’ from the spatial topographies._ + +In the previous section it was briefly discussed that the spatial fingerprint +of a given (dipolar) source, its leadfield, can be uncoupled from its moment, +i.e., the amplitude and orientation. In other words, we can construct for a +given location a leadfield matrix (consisting of 3 columns) that defines the +spatial topography of a unit amplitude source in the x, y, and z directions. +Consequently, we can compute any spatial topography due to a dipole at that +location with an arbitrary strength and orientation by right multiplying this +leadfield matrix by an appropriate dipole moment vector: `x = L*m`, where +m consists of three elements. + +To take this even further, rather than defining m as a 3x1-vector, m can be +defined as a 3xt matrix, where each column represents the dipole moment at a +particular time point. If the orientation of the dipole does not change over +time, the dipole is said to have a fixed orientation, and the matrix that represents +the dipole moment as a function of time can be decomposed into a 3x1-vector +representing the orientation (typically normalized to unit length) and a 1xt-vector +describing the activation time course of the dipole. In a formula: `y = m*s`. +If the orientation of the dipole changes over time, e.g. when the activation +is spreading over the crown of a gyrus into the bank of a sulcus, the dipole +is said to be a rotating dipole. In that case, the latter simplification is +not appropriate. + +Let’s explore these issues in the next exercises. We begin with the construction +of some sensor data that measures the activity of a single, fixed orientation +dipole. To this end, we first need to compute a leadfield: + + dippos = [0 0 8]; + sensors = ni2_sensors('type', 'eeg'); + headmodel = ni2_headmodel('type', 'spherical', 'nshell', 3); + leadfield = ni2_leadfield(sensors, headmodel, dippos); + +Note, that by specifying only 3 dipole parameters, we get a leadfield matrix, +consisting of 3 columns. Now, we create an activation time course: + + [sourceactivity, time] = ni2_activation('ncycle', 0.5, 'frequency', 1); + +We now apply an orientation vector dipmom to the activation time course, and +then apply the leadfield to this. All this can be done in a single step: + + dipmom = [1; 0; 0]; + sensordata = leadfield * dipmom * sourceactivity; + +**_Q4.1 - Take some time to understand this multiplication. Also, verify that +this yields the same results as including the dipole moment parameters in the +call to `ni2_leadfield`._** + +These simulated sensor data are actually not interesting at all. This can +be seen when visualizing the sensor activations in a movie: + + ni2_topomovie(sensors, sensordata, time); + +**_Q4.2 - Describe how the topography changes over time. Make a distinction +between the ‘shape’ of the topography, and the ‘strength’._** + +You have hopefully observed that the shape of the topography does not change +over time. If we collapse the leadfield with the dipmom in the first equation +on this page with + + leadfield1 = leadfield * dipmom; + whos leadfield leadfield1 dipmom + +the first equation can be written as: + + sensordata1 = leadfield1 * sourceactivity; + whos sensordata1 + +Where both `leadfield` and `sourceactivity` are vectors (a column and a row +vector respectively). Consequently, each row in `sensordata` can be considered +a scaled version of the vector `sourceactivity`. Also, each column in `sensordata` +can be considered a scaled version of the vector `leadfield`. Thus, the `sensordata` +matrix is a bit boring. + +Once the spatial topographies consist of more than one column (and consequently +there are multiple source activation time courses), the structure in the spatio-temporal +data matrix becomes much richer. For example, we can construct a "rotating" +dipole in the following way: + + sourceactivity2 = ni2_activation('ncycle', 0.5, 'frequency', 1, 'latency', 0.6); + sourceactivity3 = zeros(size(sourceactivity2)); + + sourcedata = [sourceactivity; sourceactivity2; sourceactivity3]; + figure; plot(time, sourcedata); legend({'1', '2', '3'}) + + sensordata2 = leadfield * sourcedata; + ni2_topomovie(sensors, sensordata2, time); + +**_Q4.3 - Describe how the topography changes over time._** + +In the previous exercise we have a constructed a rotating dipole from 2 activation +time courses, where the leadfields corresponding to each of the time courses +were defined at a single location (with different orientation). + +From this situation it is just a small step to consider neural activations +that are originating from different locations. The underlying mathematics is +exactly the same. + + sourceactivity3 = ni2_activation('ncycle', 0.4, 'frequency', 1, 'latency', 0.6); + sourcedata = [sourceactivity; sourceactivity2; sourceactivity3]; + whos sourcedata + + dippar = [0 0 8.5 0 1 0; 3 5 6 0.7 0 -0.7; -5 6 3 0.5 -0.5 0.6]; + leadfield = ni2_leadfield(sensors, headmodel, dippar); + whos dippar leadfield + sensordata3 = leadfield * sourcedata; + whos leadfield sourcedata sensordata3 + ni2_topomovie(sensors, sensordata3, time); + +**_Q4.4 - Describe the changes in topography over time._** + +**_Q4.5 - Create sensor data using three or four source with different activation +patterns to get a feel for the mixing process. In the next class we will compare +the activation patterns that you have generated, and take a vote on the most +aesthetic spatiotemporal movie._** + +## 5 Differences between EEG and MEG, and the issue of the EEG reference + +This section deals with the fundamental differences between EEG and MEG, which +is threefold: # EEG picks up more activity than MEG, # EEG signals are more +blurred than MEG, and # with EEG one always has to consider the reference electrode. + +After doing these exercises: + +- You understand the difference between EEG and MEG in the spatial topography for a given source. +- You understand that MEG does not see all activity +- You understand the influence of the reference electrode on the spatiotemporal activation pattern. + +We start with comparing spatial topographies generated by the same dipole, +but measured either with EEG or with MEG. + +As usual we need the following ingredients: a description of a set of sensors, +a volume conduction model, and a set of source parameters. First, we do the +EEG case: + + dippos = [0 0 6]; + dipmom = [1 0 0]; + dippar = [dippos dipmom]; + + headmodeleeg = ni2_headmodel('type', 'spherical', 'nshell', 3); + sensorseeg = ni2_sensors('type', 'eeg'); + leadfieldeeg = ni2_leadfield(sensorseeg, headmodeleeg, dippar); + +We can plot the topography in 2D using: + + ni2_topoplot(sensorseeg, leadfieldeeg); + +or plot everything in 3D with: + + figure + ft_plot_headmodel(headmodeleeg, 'facealpha', 0.3) + ft_plot_sens(sensorseeg, 'elecshape', 'circle', 'label', 'on') + ft_plot_dipole(dippos, dipmom) + ft_plot_axes([], 'unit', 'cm') + ft_plot_topo3d(sensorseeg.chanpos, leadfieldeeg) + ft_colormap('*RdBu'); + title('EEG') + +Then the MEG case: + + headmodelmeg = ni2_headmodel('type', 'spherical', 'nshell', 1); + sensorsmeg = ni2_sensors('type', 'meg'); + leadfieldmeg = ni2_leadfield(sensorsmeg, headmodelmeg, dippar); + +Again, the topography can be plotted with: + + ni2_topoplot(sensorsmeg, leadfieldmeg); + +and in 3D with: + + figure + ft_plot_headmodel(headmodelmeg, 'facealpha', 0.3) + ft_plot_sens(sensorsmeg, 'coilshape', 'circle', 'label', 'on') + ft_plot_dipole(dippos, dipmom) + ft_plot_axes([], 'unit', 'cm') + ft_plot_topo3d(sensorsmeg.chanpos, leadfieldmeg) + ft_colormap('*RdBu') + title('MEG') + +Note that for the MEG situation, the volume conduction model can consist of +just a single sphere and the sensors are just above the head surface, not directly +touching. + +**_Q5.1 - Describe and explain the differences in topography._** + +**_Q5.2 - Explore the effect of the source parameters on the MEG topographies. +To this end, you can vary the location of the dipole, as well as the dipole +moment. In particular:_** + +1. **_Observe what happens if you change the depth of the dipole, i.e., by moving +the position closer to the center of the head._** +1. **_Observe what happens if the dipole moment has a non-zero radial component, +where in the extreme case the orientation is fully radial. See below for some +notes with respect to dipole orientation._** + +**_Q5.3 - Explore what happens if there is more than one source in the topography. +For example, model 2 dipoles with positions [4 0 6] and [-4 0 6] both with a +moment [1 0 0], and visualize the combined leadfield. Use matrix multiplication +to do the instantaneous mixing._** + +In a spherical volume conductor model, a dipole that has its orientation aligned +with the radius of the sphere is called a radial dipole. More generally, if +you draw a line between the origin of the sphere and the position of the dipole, +the projection of the dipole moment onto this line is called the radial component +(see figure). The plane perpendicular to the radius is called the tangential +plane. An MEG-sensor is blind to the radial component of any dipole in a spherical +volume conductor. In the most extreme case, the dipole moment is fully radially +oriented, yielding to a zero magnetic field. In reality, most dipoles have both +a non-zero radial component, as well as 2 tangential components. Also, realistically +head-shaped volume conductors are not spherical, hence a true 'radial’ orientation +does not exist. + +One very important issue in EEG data analysis is the notion that we are dealing +with potential differences. Typically, however, the signals are thought of as +representing the electric potential measured at a particular location on the +scalp. In reality the signals represent the difference in electrical potential +measured at these particular locations, and the potential measured at a (distant) +reference electrode. In an ideal situation the reference electrode is electrically +silent, but this is never true. This is due to intrinsic electrode noise, and +due to the fact that the reference electrode will also pick up some biological +signal. In most practical situations, a reference electrode at the location +of one of the mastoids is used. This location for sure picks up potential fluctuations +due to brain activity. Therefore, EEG measurements are rereferenced to a linked-mastoid +reference. As a consequence each signal represents the difference between the +potential at each location and the average potential at the mastoid electrodes. +This is not really an optimal referencing scheme, but it is widely used and +thus facilitates comparison of results with previous literature. + +We will now explore the effect of the reference on the spatiotemporal activations +measured at the scalp. First we will simulate some EEG data using a single source +that has some activation time course. By now the following code should start +to look familiar to you. + + headmodel = ni2_headmodel('type', 'spherical', 'nshell', 3); + sensors = ni2_sensors('type', 'eeg'); + leadfield = ni2_leadfield(sensors, headmodel, [0 0 6 1 0 0]); + + [sourceactivity, time] = ni2_activation; + sensordata = leadfield * sourceactivity; + +Note the last line of code above. Remember that this represents the application +of the formula for the basic observation model. + +The algorithm that computes the leadfield also has to assume a reference electrode. +In this case, an average reference is assumed, meaning that each of the values +in the data represents the potential difference between the potential at the +channel’s location and the average potential across all channels. + +**_Q 5.4 - Confirm that the average over all channels is zero by plotting the +average versus the time. You might observe that it is slightly different from +zero; that is because computers have a limited numerical precision._** + +In order to emulate the situation that occurs just when the data has been +recorded, we need to rereference these data. As the reference we pick a channel +at the back of the head, just posterior of where the negativity has its maximum. + +As before, we can visualize the topography of the activity at any time point, +using the plot function. Let’s use the 485th timepoint here and let's add some +options in the "cfg" argument to show the channel numbers: + + cfg = []; + cfg.marker = 'numbers'; + + ni2_topoplot(sensors, sensordata(:, 485), cfg); colorbar + +Re-referencing is nothing else than subtracting two signals. This can be easily +done in MATLAB: + +channel 74 is an occipital channel + sensordata_reref = sensordata - repmat(sensordata(74,:), [91 1]); + + ni2_topoplot(sensors, sensordata_reref(:, 485)); colorbar + +**_Q5.5 - Also visualize the rereferenced data at the same time point and +describe the similarities/differences. Hint: don’t forget to look at the values +that are represented by the colors, using the colorbar._** + +Now have a look at the time courses, by doing: + + figure; plot(time, sensordata); + figure; plot(time, sensordata_reref); + +**_Q5.6 - Describe the differences._** + +A common rereferencing strategy is to use the linked mastoids as a reference. +This is achieved by taking the data (remember: the data has been physically +referenced during the recording to one of the mastoid electrodes), and subtract +half of the signal recorded at the other mastoid electrode. In MATLAB it looks +like this: + + sensordata_reref2 = sensordata - repmat(mean(sensordata([69 79],:),1), [91 1]); + +Note that we take the data as referenced to the 'left mastoid’ because this +represents more realistically the data how it would be recorded. Note also, +that our 'right mastoid’ corresponds with channel number 87. + +**_Q5.7 - Visualize the rereferenced to linked mastoids data and compare it +with the other results._** + +**_Q5.8 - Repeat the steps above, but now with a dipole that is close to the +‘left mastoid’ electrode, e.g. using the dipole parameters `[0 7 1 0 1 0]`._** diff --git a/workshop/neuroimaging2-2526/sr3.md b/workshop/neuroimaging2-2526/sr3.md new file mode 100644 index 000000000..bd75761e4 --- /dev/null +++ b/workshop/neuroimaging2-2526/sr3.md @@ -0,0 +1,519 @@ +--- +title: SR3 - Source reconstruction +tags: [neuroimaging2-2425] +--- + +## 6 Fitting EEG/MEG activity with dipole models + +This document contains the MATLAB exercises that form part of the course “Neuroimaging +II” relating to dipole modeling, or dipole fitting. + +The best way to approach the exercises is to create a MATLAB-script that can +be used to copy-and-paste the MATLAB-code in the instructions, and allows you +to play around with various parameter settings etc. You can open a new script +in the MATLAB-editor by typing edit on the command line, or by going to File->New->Script. +Next, you can save each script as a separate file in your course folder. + +Set up the MATLAB path + + ni2_startup + ft_version + +_These exercises treat the topic of dipole modeling. This technique is a so-called +parametric inverse method, because we will estimate a set of parameters yielding +a model of the observed data, which is aimed at explaining the data in an optimal +way. The modeling assumption we make here is that the spatial distribution of +measured potential or magnetic field at any moment in time is 'caused’ by a +limited number of equivalent current dipoles (ECDs). In the past sessions you +have learnt that we can build models of the potential difference or magnetic +field distribution for a given set of dipole parameters (3 location, and 3 moment +parameters). This is the so-called forward model. The idea of dipole modeling +is as follows: given a set of parameters we build a forward model of the data, +and compare this forward model with the measured data, by computing a measure +of goodness-of-fit. Then, we manipulate the parameters a tiny bit, recompute +the forward model, and see whether the goodness-of-fit is improved. If the goodness-of-fit +is improved, then the current model of the data apparently is better than the +previous one. We continue changing the parameters until we cannot further improve +the goodness-of-fit._ + +### 6.1 Fitting a model to the observed data + +The first thing we need to do is to determine which data matrix we are going +to model. As you know by now, the measured data can be represented in a spatio-temporal +matrix, where the rows represent the channels, and the columns represent the +spatial topography of the observed activity for a given time point. We will +start by modeling the data that is observed at a single time point. Then, we +need to make an assumption about the number of dipoles (and thus the number +of parameters) that we believe underly the observed data. Here, we start with +assuming that the data can be modeled with just a single dipole. + +After doing these exercises + +- You understand how the observed data can be fitted to a model. +- You understand how to quantify the goodness-of-fit + +We begin by generating some sensor-level data, by now this hopefully starts +to look familiar: + + [data, time] = ni2_activation; + sensors = ni2_sensors('type', 'eeg'); + headmodel = ni2_headmodel('type', 'spherical', 'nshell', 3); + leadfield = ni2_leadfield(sensors, headmodel, [4.9 0 6.2 0 1 0]); + sensordata = leadfield * data + randn(91, 1000)*1e-3; + +We have added some noise to the data, in order to make it look more realistic. + + figure; plot(time, sensordata); + +Now we select a single time slice from the data matrix, which will serve as +our observed topography: + + topo_observed = sensordata(:, 500); + figure; ni2_topoplot(sensors, topo_observed); + +**_Q6.1.1 - Describe where in the scalp topography the dipolar pattern is visible._** + +The standard-error-of-mean (SEM) or noise in the average scales inversely +proportional to the square root of the number of trials. So with 10 times more +trials, the noise is a factor sqrt(10) smaller, which is about 3x. + +**_Q6.1.2 - Repeat the process with 3x times as little noise and make a figure. +This corresponds to 10x more trials._** + +**_Q6.1.3 - Repeat the process with 3x times as much noise and make a figure. +This corresponds to 10x fewer trials._** + +What is typically done in dipole fitting is that the model of the observed +data is created in a two-step procedure. First, given a set of location parameters, +the forward model is created, yielding an Nchannels x (Ndipoles x 3) leadfield +matrix, where each triplet of columns in the leadfield matrix represents the +topography of a unit-amplitude dipole at a given location, with an orientation +along one of the axes of the coordinate system. In a second step, the parameters +that describe the dipole(s’) moment are estimated using linear estimation. This +can be seen as a linear regression problem, where the leadfield matrix is the +'design matrix’, and the observed topography is the dependent variable. + +In MATLAB this can be done in the following way. We assume a dipole at an +arbitrary location and first model the leadfield. + + leadfield = ni2_leadfield(sensors, headmodel, [5 5 4]); + +Next, we are going to estimate the remainder of the parameters, i.e., the +3 dipole moment parameters, using linear estimation/least-squares regression. +In MATLAB this can be easily done using the `\` or backslash operator. + +since we have + + topo_opserved = leadfield * dipmom + noise + +By ignoring the noise we get + + topo_opserved = leadfield * dipmom + +or by flipping left and right side + + leadfield * dipmom = topo_opserved + +We can (conceptually) divide both sides by the leadfield to solve for dipmom. +However, since we are not working with scalars but with matrices, the division +is not simply a scalar division. + +**_Q6.1.4 - Look up the help for the / (slash) and the \ (backslash) operator +by typing "doc /" and "doc \"._** + +We can compute the least-squares optimal solution for dipmom using + + dipmom = leadfield \ topo_observed; + +The modeled data can now easily be obtained by multiplying the leadfield +with the estimated dipole moment. + + topo_modeled = leadfield * dipmom; + +**_Q6.1.5 - Plot the topography of the modeled dipole and compare it to the +previously plotted observed topography. Describe the difference._** + +Now we can quantify the difference between the observed and modeled data, +by summing the squared differences, and relating this number to the sum of the +squared data values. The smaller this number, the better the fit. + + topo_residual = topo_observed-topo_modeled; + sumsq = sum(topo_residual.^2) ./ sum(topo_observed.^2); + +The value we have obtained for sumsq does not yet mean that much, unless it +is compared to the sumsq obtained for a fitted model with different parameters, +i.e., a dipole at another location. + +**_Q6.1.6 - Make a scatter plot of the values in topo_modeled (along x) and +topo_observed (along y)._** + +**_Q6.1.7 - Compare the sumsq value to the correlation that you can observe +(and compute) between the topo_modeled and topo_observed._** + +### 6.2 Finding the optimal model + +Now we know how to model the observed data using the leadfields created for +one or more dipoles with a prespecified location, and we know how to quantify +the goodness-of-fit between the modeled and observed data. Next, we need to +consider the strategies that can be used to find the optimal model. Since the +leadfields are non-linear functions of the parameters, there is no easy analytic +solution to this problem. Therefore, the implicit strategy is to sample the +parameter space, and to quantify for each setting of the parameters the goodness-of-fit. +The parameter settings yielding the best model fit are selected. Typically, +it does not make sense to just start placing dipoles at random locations and +to hope that you will find the model with the best overall fit (out of all possible +models. Particularly, when more than one dipole is assumed, the number of potential +combinations of locations quickly explodes, and becomes unmanageable. However, +in the single dipole case this seems to be a reasonable strategy. Under the +assumption that the error landscape (i.e., 1 – goodness-of-fit expressed as +a function of dipole location) is relatively smooth, one can sample the total +set of possible source locations on a 3-dimensional grid of equally spaced dipole +locations, and select the location that yields the smallest error. From this +location, one could start a non-linear search to find a final solution. + +After these exercises + +- You will understand the concept of the error function +- You will understand the concept of a grid search + +To create an error function, one simply needs to repeat multiple times the +model fitting described in the previous section. One sensible strategy is to +create a regular 3-dimensional grid of dipole positions that can be used to +sample the parameter space. Here we do this with the following function call: + + sourcemodel = ni2_sourcemodel('type', 'grid', 'headmodel', headmodel, 'resolution', 1); + +As a little aside, this creates a variable in MATLAB that is a so-called structure, +which is a special type of variable that is convenient when working with data +objects that have multiple attributes, i.e., features that belong together. +These features are stored in so-called fields, which contain the actual data. +In this example, if you type `sourcemodel` on the command line, you will see: + + sourcemodel = + pos: [3610x3 double] + inside: [3610x1 logical] + dim: [19 19 10] + +which means that the structure called `sourcemodel` has 4 fields, called `pos`, +`inside`, `outside`, and `dim`. The most relevant field is the pos-field which +contains a matrix with the positions defined in 3D-space. Note, that the number +of positions is 3610, which equals the product of the elements in the dim-field. +The dim-field defines the number of dipoles in each of the three directions +of the 3D-box that defines the search space. The inside and outside fields are +vectors which index which of the dipole positions are within the volume conduction +model, and which ones are outside it. Note, that the total number of elements +in the inside and outside fields matches the total number of positions. FieldTrip +makes extensive use of MATLAB-structures, so it is good to have some basic understanding +of this type of variable. + +For illustration purposes we will first look at a subset of all positions +in this 3D-grid, and select those positions which have a z-coordinate of 6. +In other words, we are selecting those points that are on a plane that goes +almost through the location where we actually simulated the activity (we are +cheating a bit here, because usually we don’t know this of course). + + pos = sourcemodel.pos(sourcemodel.inside,:); + sel = find(pos(:, 3)==6); + pos = pos(sel,:); + +**_Q6.2.1 - How many source positions do we have remaining after selecting +only the positions inside the head, and only for this slice?_** + +Now, what we can do is repeat the steps in the previous section for each of +these points, i.e., we will model for each of the positions a dipole that optimally +explains the observed topography and compute a measure of goodness-of-fit. We +store these goodness-of-fit measures in a vector, so that we later can determine +which position gave the best fit. Doing the same time multiple times can be +easily solved with a for-loop. + + sumsq = zeros(size(pos, 1), 1); + + for k=1:size(pos, 1) + disp(k) + leadfield = ni2_leadfield(sensors, headmodel, pos(k,:)); + dipmom = leadfield \ topo_observed; + topo_modeled = leadfield * dipmom; + topo_residual = topo_observed-topo_modeled; + sumsq(k) = sum(topo_residual.^2) ./ sum(topo_observed.^2); + end + +Now we have a variable sumsq that is a vector, rather than a single number. +We can visualize this variable using the `plot` function. + +**_Q6.2.2 - Plot the sumsq values (along the y-axis) versus the index of the +dipole (along the x-axis)._** + +We can now look for the position that yields the lowest sum of squared difference +values, i.e., the one with the best model fit. We could achieve this by zooming +in into the figure and writing down the x-coordinate where the sumsq-variable +seems to have the lowest value, but we can also use MATLAB’s min function: + + [m, ix] = min(sumsq) + pos(ix,:) + +Depending on the noise in the data, this position will be more or less close +to the actual simulated location. + +One can improve on this estimate by fitting a model to more than one time +point. The idea is that the influence of noise on the single time point observed +topographies is attenuated when combining across time points. Assuming a fixed +dipole location, this extended temporal model is quite straightforward, we will +return to this later. + +First, for convenience, we will repeat the steps above by using a set of pre-computed +leadfields, that we compute for the *full 3D grid*. Since the computation of +the leadfields is relatively expensive (time-consuming) in terms of computation, +it makes sense to compute them only once when they will be used multiple times. +Further down in this tutorial we will explore the effect of model assumptions +on the outcome, and do this using the grid search, so it pays off to compute +the leadfields. + + if false + pos = sourcemodel.pos(sourcemodel.inside,:); + leadfield = ni2_leadfield(sensors, headmodel, pos); + else + load leadfields + end + + whos leadfield pos + +Now we can very quickly compute the goodness-of-fit for all locations in the +3D grid: + + sumsq = ones(size(pos, 1), 1); + + for k=1:size(pos, 1) + disp(k) + ik = (k-1)*3+(1:3); % select the three columns corresponding to one position + dipmom = leadfield(:, ik) \ topo_observed; + topo_modeled = leadfield(:, ik) * dipmom; + topo_residual = topo_observed-topo_modeled; + sumsq(k) = sum(topo_residual.^2) ./ sum(topo_observed(:).^2); + end + +### 6.3 Adding the time dimension + +When using the topography from a single time point for dipole modeling, noise +in the data can negatively affect the result. To this end it makes sense to +combine data across time points. Under the assumption that the model is stationary +in terms of dipole location (and possibly also in terms of orientation), pooling +across time points can be efficient in improving the model fit. + +After these exercises: + +- You will appreciate the fact that including temporal information can improve +the dipole fit. + +For this exercise we will the simulated data that was generated in section +2. If you are not working through these exercises in a single session, please +go back to section 2, and re-create the variable sensordata. + +Previously, we have used time point '500’ for the extraction of the observed +topography. + +Yet, in our very simple generative model (with generative model we mean the +model that actually underlies the simulated data) there is only one dipole with +a fixed location and orientation. In other words, and this is something that +was highlighted in a previous session, there is essentially only a single topography +present in the data, which only varies as a function of time. Thus, in principle, +using time point '499’ should yield the same result as using time point '501’, +or, for that matter, any other time point. Obviously, since the strength of +the activity varies over time, some time points are more informative than others. + +Now we are going to fit a single dipole model to a spatial topography from +a single time point, just as in the previous section. To speed things up, we +will not compute the leadfields on the fly, but load the pre-computed leadfields. + + topo_observed = sensordata(:, 500); % select one timeslice + + pos = sourcemodel.pos(sourcemodel.inside,:); + load('leadfields'); + + sumsq = ones(size(pos, 1), 1); + for k=1:size(pos, 1) + disp(k) + ik=(k-1)*3+(1:3); + dipmom = leadfield(:, ik) \ topo_observed; + topo_modeled = leadfield(:, ik) * dipmom; + topo_residual = topo_observed-topo_modeled; + sumsq(k) = sum(topo_residual(:).^2)./sum(topo_observed(:).^2); + end + [m, ix] = min(sumsq) + pos(ix,:) + +Note the use of `topo_residual(:)` in the line of code just before the end +statement. This `(:)` instruction tells MATLAB to reshape a matrix into a vector. +This is not needed in case the topography consists of just a single time point +(because the observed topography is already a vector), but we will need it when +fitting multiple time points. + +**_Q6.3.1 - Repeat the computation with topo_observed corresponding to all +sensordata, i.e., not selecting the 500th column but taking the whole matrix._** + +**_Q6.3.2 - What is the size of the topo_residual variable?_** + +**_Q6.3.3 - What is the size of the sumsq variable? Explain this._** + +### 6.4 The real deal + +So far, for didactical purposes, we have constrained ourselves and evaluated +the error-function based on a grid search. In this section we are going to use +one of FieldTrip’s core functions `ft_dipolefitting` to do a proper dipole fit +on our simulated EEG-data. + +After this section + +- You understand that the non-linear search for optimal dipole parameters +can yield a solution very close to the ground truth. + +We start again by ensuring that we have the following variables in the MATLAB-workspace: +sensordata, headmodel and sensors. If you don’t have these variables available, +please go back to section 2 and create these variables. + +Next, we need to create some variables that are recognized by FieldTrip. Most +FieldTrip functions work by providing 2 input arguments, a cfg-structure that +contains parameters that determine the behaviour of the function, and one (or +more) data-structures providing the data on which the function operates. + +Let’s first create the data variable: + + data = []; + data.avg = sensordata; + data.time = time; + data.label = sensors.label; + data.elec = sensors; + data.dimord = 'chan_time'; + +And the cfg-structure: + + cfg = []; + cfg.gridsearch = 'yes'; + cfg.model = 'regional'; + cfg.headmodel = headmodel; + cfg.latency = [0.49 0.51]; + cfg.nonlinear = 'yes'; + cfg.numdipoles = 1; + +Now we can call the function: + + dip = ft_dipolefitting(cfg, data); + +The output variable to `ft_dipolefitting` has a field `dip` containing information +about the optimal model. In particular, have a look at `dip.dip.pos`, and `dip.dip.mom`. +We can also visualize the modeled topography, and compare this to the observed +topography. These are represented in `dip.Vmodel` and `dip.Vdata`, respectively. + +**_Q6.4.1 - Using the output of ft_dipolefitting, plot a topography of the +observed data at 500ms, plot a topography of the model data at 500ms, and plot +the topopgraphy of the residual._** + +### 6.5 It’s all about the assumptions + +One important issue in dipole modeling is that the prior assumptions critically +constrain the final model (and thus the model fit). If these assumptions don’t +coincide with what’s actually in the data, this can lead to erroneous interpretations. +This can work against you in two directions, either your model is too simplistic +(i.e., you assume too few dipoles), or too complicated (you assume too many +dipoles). + +After these exercises: + +- You understand that if your model assumptions violate the underlying data, +the results are suboptimal. + +First, we create some sensor-level data that contains two sources: + + sensors = ni2_sensors('type', 'eeg'); + headmodel = ni2_headmodel('type', 'spherical', 'nshell', 3); + + [data1, time] = ni2_activation; + [data2, time] = ni2_activation('frequency', 11, 'latency', 0.48); + + leadfield1 = ni2_leadfield(sensors, headmodel, [4.9 0 6.2 0 1 0]); + leadfield2 = ni2_leadfield(sensors, headmodel, [-5.3 0 5.9 1 0 0]); + sensordata = leadfield1*data1 + leadfield2*data2; + +**_Q6.5.1 - Plot the EEG data versus time (along the x-axis) and observe the +alternating pattern._** + +**_Q6.5.2 - Use the `ni2_topomovie` function to explore the alternating pattern +of activity._** + +Let's add some noise + + sensordata = leadfield1*data1 + leadfield2*data2 + randn(91, 1000)*.7e-3; + +Now we will perform the same grid search as in section 4, using the 490th +time point for the observed topography: + + topo_observed = sensordata(:, 490); + + sourcemodel = ni2_sourcemodel('type', 'grid', 'resolution', 1); + pos = sourcemodel.pos(sourcemodel.inside,:); + load('leadfields'); + + sumsq = ones(size(pos, 1), 1); + for k=1:size(pos, 1) + disp(k) + ik=(k-1)*3+(1:3); + dipmom = leadfield(:, ik) \ topo_observed; + topo_modeled = leadfield(:, ik)*dipmom; + sumsq(k)= sum((topo_observed(:)-topo_modeled(:)).^2) ./ sum(topo_observed(:).^2); + end + +Now we are going to use the FieldTrip `ft_dipolefitting` function to fit a +single dipole: + + data = []; + data.avg = sensordata; + data.time = time; + data.label = sensors.label; + data.elec = sensors; + data.dimord = 'chan_time'; + + cfg = []; + cfg.gridsearch = 'no'; + cfg.model = 'regional'; + cfg.headmodel = headmodel; + cfg.latency = [0.49 0.51]; + cfg.nonlinear = 'yes'; + cfg.numdipoles = 1; + dip = ft_dipolefitting(cfg, data); + +**_Q6.5.3 - Where did the fitted dipole end up and how does that compare to +the initial position of dipole 1 and dipole 2?_** + +We can also fit a model with two dipoles, this can be easily achieved by changing +the `cfg.numdipoles` option into 2. + + cfg = []; + cfg.gridsearch = 'no'; + cfg.model = 'regional'; + cfg.headmodel = headmodel; + cfg.latency = [0.49 0.51]; + cfg.nonlinear = 'yes'; + cfg.numdipoles = 2; + dip = ft_dipolefitting(cfg, data); + +As you may have noticed, the result is not particularly accurate. The reason +for this is that the optimization algorithm got trapped in a local minimum of +the error function. This is more likely to happen, the more complicated the +underlying model (i.e., more free parameters lead to a high-dimensional error +function with a complicated structure and potentially many local minima). We +can however inform the fitting algorithm with dipole positions from which to +start the non-linear search. If these starting positions are sufficiently close +to the actual source positions, the algorithm will converge to the correct solution. + + cfg = []; + cfg.gridsearch = 'no'; + cfg.model = 'regional'; + cfg.headmodel = headmodel; + cfg.latency = [0.49 0.51]; + cfg.nonlinear = 'yes'; + cfg.numdipoles = 2; + cfg.dip.pos = [5 0 5; -5 0 5]; % starting positions + dip = ft_dipolefitting(cfg, data); + +**_Q6.5.4 - Repeat the dipole fit with two dipoles with their initial positions +for the whole tine interval from 0 to 1000 ms. Plot the residual variance (dip.dip.rv) +versus time and describe what you see._** diff --git a/workshop/neuroimaging2-2526/sr4.md b/workshop/neuroimaging2-2526/sr4.md new file mode 100644 index 000000000..bddc371ff --- /dev/null +++ b/workshop/neuroimaging2-2526/sr4.md @@ -0,0 +1,544 @@ +--- +title: SR4 - Source reconstruction +tags: [neuroimaging2-2425] +--- + +## 7 Modelling EEG/MEG activity using distributed sources + +This document contains the MATLAB exercises that form part of the course “Neuroimaging +II” relating to the minimum norm inverse methods for underdetermined systems. + +The best way to approach the exercises is to create a MATLAB-script that can +be used to copy-and-paste the MATLAB-code in the instructions, and allows you +to play around with various parameter settings etc. You can open a new script +in the MATLAB-editor by typing edit on the command line, or by going to File->New->Script. +Next, you can save each script as a separate file in your course folder. + +Set up the MATLAB path + + ni2_startup + ft_version + +_If we have measured data and a pre-computed leadfield (based on known sensor +positions and a volume conductor model), then we would like to compute an estimate +of the sources. Many methods use a linear matrix to compute this 'inverse solution’, +as a linear weighting of the senors to compute activity at a source location. +Throughout this homework, the equation representing the data and model is `y=L*s+N`. +Here, `y` is the sensor data, `L` is the known leadfield, `s` is the source +amplitude or dipole moment and N is the sensor noise._ + +_With distributed source models and the minimum norm solution the observed +data is assumed to be determined by a mixing of many active sources, where the +number of active sources is much higher than the number of observed channels. +This leads to an underdetermined problem, and one way to solve this is to assume +that a good guess of the distribution of the source amplitudes has the lowest +possible norm. One important feature of distributed source models is the fact +that the amplitudes of the individual sources are estimated in a single inversion +step._ + +### 7.1 Underdetermined linear systems of equations + +For non-unique solutions (underdetermined systems) we can use the minimum +norm of the overall source strength as additional constraint. + +After this section, you will + +- Understand that multiple solutions exist for an underdetermined system +- Understand that the pseudo-inverse of the leadfield matrix gives the solution +with the minimum-norm of the source power + +First begin with a 'toy’ composite leadfield matrix. It represents 3 source +locations and 2 sensors. + + L = [1 0 1; + 0 1 1]; + +Although it is an unrealistic example, it is useful to gain some intuition. +Let’s assume the data `y` we measured at a single time point is + + y = [2 1]'; + +Now we want to solve for what the source amplitude `s` is, given this data +`y` and given the known leadfield `L`. + +Assume no noise `N` for now, so the equation is `y = L*s`. + +This matrix multiplication essentially is a linear system of equations, where +the number of equations is 2 (in general: the number of rows in the leadfield +matrix), and the number of unknowns is 3 (in general: the number of columns +in the leadfield matrix. This never has an unique solution; there are many solutions +that satisfy this equation. Let’s first explicitly rewrite the matrix multiplication. +Refresh your understanding of matrix multiplication if needed. + + y(1) = L(1, 1)*s(1) + L(1, 2)*s(2) + L(1, 3)*s(3); % equation 1 + y(2) = L(2, 1)*s(1) + L(2, 2)*s(2) + L(2, 3)*s(3); % equation 2 + +Filling in the numbers, we get: + + y(1) = 1*s(1) + 0*s(2) + 1*s(3) = 2 + y(2) = 0*s(1) + 1*s(2) + 1*s(3) = 1 + +where s(1), s(2) and s(3) are the unknowns. + +Note that equation 1 represents a matrix multiplication of the first row of +the leadfield matrix with the source vector, and that equation 2 represents +the matrix multiplication of the second row of the leadfield matrix with the +source vector. Note, also, that we have more unknowns (3) than equations (2). +In other words, we could take an arbitrary value for, say, s(3), and we will +still be able to find a valid solution to the linear system of equations. + +By analogy, you may remember from high school geometry, that equations with +3 unknowns describe a plane in 3-dimensional space, and that 2 planes (if they +are not parallel) intersect in a line. This line represents all valid solutions +to the linear system of equations. + +In this toy example, it is very straightforward to parameterize the solution. +Let’s call the value that we take for s(3) the value `a`. Using this substitution, +and applying it to the equations above, we get: + +the variable a can be anything, so let's give it a random value + a = rand(1) + +note that vector s should be a column-vector, not a row-vector + s(1,1) = 2-a; + s(2,1) = 1-a; + s(3,1) = a; + +**_Q7.1.1 - Show that this solution for `s` satisfies the linear system of +equations when a=0._** + +**_Q7.1.1 - Show that this solution for `s` satisfies the linear system of +equations when a=100._** + +We will use this parameterization in the next section. + +### 7.2 The 'best’ solution based on additional constraints: minimize norm. + +In the example presented above, mathematically it is equally valid to take +a value of 1 for a, as it is to take a value of 100. The former will yield source +amplitudes of (1, 0, 1), whereas the latter will yield source amplitudes of +(-98,-99, 100). + +It may be realistic to assume that the source activity measured at any given +instant by MEG/EEG is primarily due to active sources that are moderately active +(a=1). The alternative would be that all sources would be highly active (a=100). +A mathematically convenient and biologically plausible constraint is to minimize +the total power of the sources across the brain (rather than taking the amplitude, +we minimize for the amplitude squared). The sum of the power of the sources +across the brain is also called the L2-norm; the sum of the amplitudes of the +sources across the brain is represented by the so-called L1-norm. + +To compute the L2-norm, we take the square root of the sum of the squares +of each element in the vector. Here’s an example of source amplitude over 3 +source locations: + + s = [0 -1 2]'; + +The source power (L2 norm) can be computed in a variety of ways in MATLAB: + + snorm1 = sqrt(s(1)^2+s(2)^2+s(3)^2) + snorm2 = sqrt(sum(s.^2)) + snorm3 = sqrt(s'*s) + snorm4 = norm(s) + +**_Q7.2.1 - Compute and write down the L2-norm_** + +Rather than trying out random values for `a` to find the solution that yields +the lowest source norm, we can use the parameterization that we defined in the +previous section: + + snorm = sqrt(s(1)^2+s(2)^2+s(3)^2) + +Filling in the parameterization we obtained in section 2.1, we get: + + snorm = sqrt((2-a)^2+(1-a)^2+a^2); + +Rearranging the terms, we get: + + snorm = sqrt((4-4*a+a^2)+(1-2*a+a^2)+a^2); + snorm = sqrt(3*a^2-6*a+5); + +**_Q7.2.2 - Compute and write down the L2-norm for a=1 and for a=100._** + +Ignoring the square root, we need to minimize `(3*a^2-6*a+5)` which can be +done in a variety of ways, e.g. using the ABC-formula (https://en.wikipedia.org/wiki/Quadratic_formula) +to find the minimum of a parabolic equation, or by taking the derivative of +this equation and finding the point along the x-axis where the derivative of +the parabolic function is 0. + +**_Q7.2.3 - Take the function `f(a)=3*a^2-6*a+5` and compute its values for +the variable "a" ranging from -100 to +100 in small steps. Make a plot of f(a) +versus a._** + +**_Q7.2.4 - What is the smallest value of f(a), and for which value of a is +that value observed?_** + +### 7.3 Inverting leadfield matrices + +For situations with more than just a few sensors and a few sources, it is +too tedious to solve the system of equations for the unknown values by hand. +Fortunately, we can use matrix algebra to solve large systems of equations, +and let the computer do the work. + +First, by mean of a detour, let’s have a look at the following equation: + + y = 3*s + +If we want to solve this (very simple system of linear equations) for `s` +(the source amplitude), assuming that we know `y`, we simply divide each side +of the equation by 3, so we get: + + 1/3 * y = s + +Instead of dividing each side of the equation by 3, we can also say that we +are multiplying each side of the equation by 1/3. 1/3 is also known as the 'inverse’ +of 3, which can also be written as 3^-1. When we are dealing with matrices, +the same logic applies: + + Y = L * S + L^-1 * Y = L^-1 L * S + L^-1 * Y = S + +Here `L^-1*L = I`, where the matrix `I` is the identity matrix, all zeros +off the main diagonal and all ones on the diagonal; this is the matrix equivalent +of the number "1". + +The above equation however only works in a very limited number of cases. The +reason is that a matrix inverse is only defined for square matrices (i.e., in +our case we would need an equal number of channels and sources), where moreover +the leadfield matrix must fulfill the mathematical property that it is actually +invertible. + +We can use the so-called pseudo-inverse of a matrix to get a solution to the +linear system of equations. The term 'pseudo’ refers to the fact that the matrix +you get after taking the pseudo-inverse behaves a bit like an inverse matrix, +but not exactly. + +The pseudo-inverse in this case can be computed in math formulation as: + + Lpinv = L' * (L*L')^-1 + +The MATLAB pinv function achieves the same: + + help pinv + Lpinv = pinv(L) + +In contrast to a proper inverse matrix, where both `A*A^-1 = I` and `A^-1*A += I`, the pseudoinverse exists in only one direction. This can be easily seen +when inspecting `Lpinv*L` and `L*Lpinv`. + +**_Q7.3.1 - Compute and write down both `Lpinv*L` and `L*Lpinv`._** + +The first one does not result in an identity matrix. The second equation results +in an identity, because it is equivalent to `(L*L')^-1 * (L*L')`. The matrix +products between the brackets are the same, both form a square matrices and +they happen to be in general invertible, thus the product between the inverse +of that product with itself will yield `I`. + +**_Q7.3.2 - What happens if you try to compute the normal inverse (using the +inv function) of matrix L?_** + +### 7.4 Pseudo-inverse of the leadfield gives minimum norm + +The pseudo-inverse is not only a clever mathematical way of computing a new +matrix that, when multiplied with the original in the correct order, gives back +the identity matrix. It happens to be a solution to the underdetermined linear +system of equations that yields the minimum norm. + +The proof for this can be derived by first showing that it is a particular +solution; subsequently we show that there is no other solution that has a norm +smaller than the one derived using the pseudoinverse. + +We know that there are many possible solutions for the equation `y = L*s`. +Let us use a capital `S` from now on to facilitate formatting; we will use `Sx` +to indicate S with the subscript `x`. If we just take two possible solutions, +`Sm` and `Sx`, where `Sm` happens to be `L'*(L*L')^-1 * y` and `Sx` being any +random other solution, we know that both `L*Sm` and `L*Sx` yield `y`. + +Thus, `L*Sm = L*Sx`, or equivalently `L*(Sm-Sx) = 0`. + +Remember from section 2.2 that the norm of the solution can be computed from +`s'*s`. This latter equation computes for each of the sources the square of +the amplitude, and sums across sources. To actually get the norm, we would also +have to take the square root of the result, but let’s not do this for now, and +look at the squared norm. Likewise, we can do `(Sx-Sm)'*Sm`, where we compute +for each of the sources the product between the amplitude modeled in `Sm` and +the amplitude difference between `Sx` and `Sm`, and sum this across sources. +Substituting `L'*(L*L')^-1 * y` for the second `Sm` in the equation, we get: + + (Sx-Sm)' * Sm = + (Sx-Sm)' * L' * (L*L')^-1 * y = + ( (Sx-Sm)' * L') * (L*L')^-1 * y = + ( L * (Sx-Sm) )' * (L*L')^-1 * y = + 0 * (L*L')^-1 * y = 0 + +Shuffling the brackets, we focus on the first 2 terms, where we can use the +general matrix property that `A'B'=(BA)`. + +Above we already concluded that `L*(Sx-Sm)` is 0, so we can also conclude +that `(Sx-Sm)'*Sm` is 0. This is because we can fill in a 0 in the last equation, +and 0 times something else will be 0. + +Now considering the squared norm of `Sx`, which is `Sx'*Sx`, we can apply +a little trick: instead of using `Sx`, we use `((Sx-Sm)+Sm)`. The latter is +of course exactly the same as `Sx`. + + Sx' * Sx = + ((Sx-Sm)+Sm)' * ((Sx-Sm)+Sm) = + (Sx-Sm)'*(Sx-Sm) + (Sx-Sm)'*Sm + Sm'*(Sx-Sm) + Sm'*Sm + +Using our previous result `(Sx-Sm)'*Sm=0`, we get + + (Sx-Sm)'*(Sx-Sm) + (s-Sm)'*Sm + Sm'*(Sx-Sm) + Sm'*Sm = + (Sx-Sm)'*(Sx-Sm) + 0 + 0 + Sm'*Sm = + (Sx-Sm)'*(Sx-Sm) + Sm'*Sm + +or + + Sx'*Sx = Sm'*Sm + (Sx-Sm)'*(Sx-Sm) + +The last equation tells us that the squared norm of `Sx` is always the sum +of the squared norm of `Sm` _plus_ the squared norm of the difference between +`Sx` and `Sm`. + +Since squared numbers are greater than or equal to zero, we can conclude that +the norm of `Sx` is always larger than or equal to to the norm of `Sm`. Hence, +`Sm` must represent the solution with the minimum norm. + +### 7.5 The real deal - start simple with simulated data without noise + +After this section you will + +- Have a basic understanding of how the minimum norm reconstruction works in practice. +- Understand why the minimum norm reconstruction tends to overemphasize activity from superficial sources. +- Know how to counteract the tendency for overemphasizing the superficial sources. +- Understand that noise in the data projects onto the estimated sources. +- Know how to counteract the contamination of the estimated source activity by the noise. + +We start by simulating some MEG data that contains two active sources. + + [activity1, time1] = ni2_activation; + [activity2, time2] = ni2_activation('frequency', 11, 'latency', 0.48); + + sensors = ni2_sensors('type', 'meg'); + headmodel = ni2_headmodel('type', 'spherical', 'nshell', 1); + + leadfield1 = ni2_leadfield(sensors, headmodel, [ 4.9 0 6.2 0 1 0]); % close to position 2352 in grid + leadfield2 = ni2_leadfield(sensors, headmodel, [-5.3 0 5.9 1 0 0]); + + sensordata = leadfield1*activity1 + leadfield2*activity2; + +Try and understand the steps above. Pay particular attention to the parameters +of the simulated dipoles. + +We now proceed to generate a MATLAB data-structure that FieldTrip understands. +This "structure" is a collection of MATLAB-variables, organized in so-called +fields, that belong together. An important aspect of these FieldTrip data structures +is that the numeric data that is represented (in our case in the 'avg’ field) +is accompanied by all information necessary to interpret this numeric data. +For example, there is a field called 'time’, that indicates each time sample +in seconds (i.e., it maps the columns of the 'avg’ field on a physical time +axis). The 'label’ field specifies the name of each channel (and tells us which +row in the 'avg’ field belongs to which channel). + + data = []; + data.avg = sensordata; + data.time = time1; + data.label = sensors.label; + data.grad = sensors; + data.cov = eye(numel(sensors.label)); + data.dimord = 'chan_time'; + +Next we will make a source reconstruction using the 'mne’ method of FieldTrip’s +ft_sourceanalysis function. Before we can do this, we need to define our source +model, i.e., the set of locations that we assume to be active. For now we assume +that the active dipoles are distributed on a regular 3D grid, with a spacing +of 1 cm between the dipoles: + + sourcemodel = ni2_sourcemodel('type', 'grid', 'resolution', 1); + + cfg = []; + cfg.sourcemodel = sourcemodel; + cfg.headmodel = headmodel; + cfg.method = 'mne'; + cfg.mne.prewhiten = 'yes'; + cfg.mne.scalesourcecov = 'yes'; + cfg.mne.lambda = 0; + cfg.keepleadfield = 'yes'; + source = ft_sourceanalysis(cfg, data); + +Let’s now have a look at the reconstructed source activity. + +**_Q7.5.1 - What is the position of grid point 2352?_** + +For each dipole location in the distributed source model, the estimated activity +is represented in the source.avg.mom field. We can easily use the MATLAB plot +command to visualize this: + + figure; plot(source.time, source.avg.mom{2352}); legend({'Sx' 'Sy' 'Sz'}); + +**_Q7.5.2 - In which direction is the dipole moment of this source at grid +location 2352 the largest? Is that consistent with the model that generated +the data?_** + +**_Q7.5.3 - What is the index of the grid point that is the closest to the +other dipole that we used in the model to generate the data?_** + +Ignoring the dipole orientation and time course of the activity, we can plot +the spatial distribution of the dipole strength that is represented in `source.avg.pow`, +which represents the squared amplitude for each time point. + + cfg = []; + cfg.funparameter = 'pow'; + cfg.location = sourcemodel.pos(2352,:); + cfg.funcolorlim = 'maxabs'; + figure; ft_sourceplot(cfg, source); + +The initial time point of which the spatial topography is plotted is t=0, +and at that moment there is no activity at all, hence the completely green distribution. +Click on the source time course in the lower right to select a time point at +which the activity in grid location 2352 peaks. + +**_Q7.5.4 - What happens with the spatial distribution of the source power +in-between two peaks at location 2352? Explain._** + +This simple noise-less example illustrates two important things. First, activity +is 'smeared’ out over various dipole locations and orientations. Second, the +estimated activity at a location closer to the sensors than the location at +which activity was simulated has a higher amplitude than the activity estimated +at the location where activity was simulated. We will return to this feature +of the minimum norm reconstruction in a later section. + +### 7.6 Simulated data with noise + +Let’s now simulate MEG sensor data with added noise: + + [activity1, time1] = ni2_activation; + [activity2, time2] = ni2_activation('frequency', 11, 'latency', 0.48); + + sensors = ni2_sensors('type', 'meg'); + headmodel = ni2_headmodel('type', 'spherical', 'nshell', 1); + + leadfield1 = ni2_leadfield(sensors, headmodel, [ 4.9 0 6.2 0 1 0]); % close to position 2352 in grid + leadfield2 = ni2_leadfield(sensors, headmodel, [-5.3 0 5.9 1 0 0]); % close to position 2342 in grid + + sensordata = leadfield1*activity1 + leadfield2*activity2 + randn(301, 1000)*.7e-10; + +Create a FieldTrip data structure: + + data = []; + data.avg = sensordata; + data.time = time1; + data.label = sensors.label; + data.grad = sensors; + data.cov = cov(randn(301, 1000)'*.7e-10); + data.dimord = 'chan_time'; + +In the field 'cov’ we create a covariance matrix that was designed to represent +the covariance of the noise in the data. This will be a relevant item when we +will discuss noise regularisation. + + sourcemodel = ni2_sourcemodel('type', 'grid', 'resolution', 1); + +Do the source reconstruction: + + cfg = []; + cfg.sourcemodel = sourcemodel; + cfg.headmodel = headmodel; + cfg.method = 'mne'; + cfg.mne.prewhiten = 'yes'; + cfg.mne.scalesourcecov = 'yes'; + cfg.mne.lambda = 0; + cfg.keepleadfield = 'yes'; + source_noise = ft_sourceanalysis(cfg, data); + +The `cfg.mne.lambda` option was set to 0. This means that the inverse solution +fits all data perfectly, where the data not only includes the activity from +the sources of interest, but also contains noise. + +We can again make a plot of the spatial distribution of power in the brain +for each time point. + + cfg = []; + cfg.funparameter ='pow'; + cfg.location = sourcemodel.pos(2352,:); + cfg.funcolorlim = [-2 2]*1e-3; % set the color limits manually + figure; ft_sourceplot(cfg, source_noise); + +**_Q7.6.1 - Select a good latency and click around in the volume; can you find +a location that shows the expected temporal pattern of activity?_** + +We can use a non-zero lambda to compute a regularized minimum norm estimate, +where this lambda parameter is used to quantify the contribution of the noise +to the measured signals. + +The larger the value for lambda, the stronger the assumed noise. In combination +with the regularization parameter, a regularized minimum-norm estimate also +requires an estimate of the noise covariance matrix. This matrix represents +the spatial structure in the noise. The noise covariance matrix can be estimated +from the data, but experimenters sometimes also use an identity matrix. The +latter strategy assumes implicitly that each channel in the data gets the same +amount of uncorrelated noise. + +Do the source reconstruction with regularisation: + + cfg = []; + cfg.sourcemodel = sourcemodel; + cfg.headmodel = headmodel; + cfg.method = 'mne'; + cfg.mne.prewhiten = 'yes'; + cfg.mne.scalesourcecov = 'yes'; + cfg.mne.lambda = 0.5; + cfg.keepleadfield = 'yes'; + source_noise_reg = ft_sourceanalysis(cfg, data); + + cfg = []; + cfg.funparameter = 'pow'; + cfg.location = sourcemodel.pos(2352,:); + cfg.funcolorlim = 'maxabs'; + figure; ft_sourceplot(cfg, source_noise_reg) + +Another way to explore the effect of regularisation is to look at the residuals +of the model. This can be obtained in the following way: + + L = cat(2, source_noise_reg.leadfield{source_noise_reg.inside}); + S = cat(1, source_noise_reg.avg.mom{source_noise_reg.inside}); + model = L*S; + residual = sensordata-model; + +**_Q7.6.2 - Show that the residual for the source_noise condition without regularization +is zero and share your figure here. Note the limited numerical precision of +the computer that might cause some numerical inaccuracy to remain._** + +### 7.7 Minimum-norm estimates 'overestimate’ the amplitude of superficial sources + +As we have seen in the previous sections, the minimum-norm estimate has a +tendency to over-estimate the amplitude of dipoles that are close to the surface. +This feature is a direct consequence of the minimum-norm constraint. In order +to explain all measured data with a source model that has the lowest possible +norm, the deep sources will be penalized more because these need to have a strong +activation in order to be picked up by the MEG sensors in the first place. + + cfg = []; + cfg.sourcemodel = sourcemodel; + cfg.headmodel = headmodel; + cfg.method = 'mne'; + cfg.mne.prewhiten = 'yes'; + cfg.mne.scalesourcecov = 'yes'; + cfg.mne.lambda = 0.5; + cfg.keepleadfield = 'yes'; + cfg.normalize = 'yes'; + cfg.normalizeparam = 1; + source_noise_lfnorm = ft_sourceanalysis(cfg, data); + + cfg = []; + cfg.funparameter = 'pow'; + cfg.location = sourcemodel.pos(2352,:); + cfg.funcolorlim = 'maxabs'; + figure; ft_sourceplot(cfg, source_noise_lfnorm) + +**_Q7.7.1 - Explain the effect that is achieved by the leadfield normalization +(cfg.normalize='yes')._** + +**_Q7.7.2 - Is the resulting spatial distribution an accurate representation +of the (in this case known) underlying source activity?_** diff --git a/workshop/neuroimaging2-2526/sr5.md b/workshop/neuroimaging2-2526/sr5.md new file mode 100644 index 000000000..32b655e49 --- /dev/null +++ b/workshop/neuroimaging2-2526/sr5.md @@ -0,0 +1,638 @@ +--- +title: SR5 - Source reconstruction +tags: [neuroimaging2-2425] +--- + +## 8 Reconstructing EEG/MEG activity using beamforming + +This document contains the MATLAB exercises that form part of the course “Neuroimaging +II” relating to the scanning methods, specifically using beamforming. + +The best way to approach the exercises is to create a MATLAB-script that can +be used to copy-and-paste the MATLAB-code in the instructions, and allows you +to play around with various parameter settings etc. You can open a new script +in the MATLAB-editor by typing edit on the command line, or by going to File->New->Script. +Next, you can save each script as a separate file in your course folder. + +Set up the MATLAB path + + ni2_startup + ft_version + +_Beamforming is the third type of strategy that can be used for source reconstruction +of EEG/MEG data that is discussed in this course._ + +_As a brief recapitulation, the *first strategy* is the dipole fitting approach, +where the underlying assumption is that the pattern in the observed data can +be explained by a limited number of neural sources that are modeled as equivalent +current dipoles. The goal is now to find the parameters for the model that optimally +explains the observed sensor-level topography. These parameters pertain to the +number, location and orientation of the dipoles, and the optimal model is determined +by minimizing the error between the model and the observed data._ + +_The *second strategy* we discussed was distributed source modeling, where +minimum norm estimation was treated in some more detail. In this approach, the +observed data is assumed to be determined by a mixing of many active sources, +where the number of active sources is much higher than the number of observed +channels. This leads to an underdetermined problem, and one way to solve this +is to assume that a good guess of the distribution of the source amplitudes +has the lowest possible norm. One important feature of distributed source models +is the fact that the amplitudes of the individual sources are estimated in a +single inversion step._ + +_The *third strategy* that will be treated here is a so-called scanning method +(we will specifically be dealing with beamforming), where the amplitude at a +predefined set of source locations is estimated iteratively, that is, for one +source location at a time. This difference with distributed source modeling +has some important consequences, which we will return to later in these exercises. +The key feature of beamforming is that no prior assumption about the number +of sources is made. Rather, for each location that we choose to evaluate, we +ask ourselves the question: 'What would the the best estimate of the dipole +moment at this location?’ It turns out, that if we assume the underlying sources +to be uncorrelated over time we can get quite good estimates of our sources’ +activation time courses using beamforming._ + +_The 'best estimate’ in the previous question is more specifically operationalized +as the estimate that best captures the activity that can originate from that +location while optimally suppressing interference from sources at other locations. +When we use so-called adaptive beamformers for this, we need not only the biophysically +constraining forward model (i.e., the leadfield matrix), we also need the sensor-level +data, or more specifically the sensor-level covariance matrix. Before we will +start playing with the beamformer and get a feel for how and why it works, the +next section will deal with the concept of the sensor-level covariance matrix._ + +### 8.1 The sensor covariance matrix: a key ingredient for the beamformer + +After this section you will + +- Understand the concept of a covariance matrix. +- Argue that the sensor covariance matrix represents a mixture of the leadfield outer products of all pairs of active sources. +- Understand that, if the sources underlying the sensor covariance matrix are uncorrelated over time, the covariance matrix reduces to a weighted sum of the leadfield outer products of only the active sources. + +In the following we will use `X` and `Y` to represent vectors rather than +matrices, to allow for formatting an element of that vector with a subscript +as `Xi`. + +The covariance between two variables (in our case time series) is defined +as: + + cov(X, Y) = ∑ (Xi-mean(X)) * (Yi-mean(Y)) / (n-1) + +If, for simplicity, we assume the vectors x and y to have a mean value of +0, this equation reduces to: + + cov(X, Y) = ∑ (Xi * Yi) / (n-1) + +In words: take the element-wise product of the elements in the vectors x and +y, and sum these products. Finish by normalizing with `(n-1)`, where n is the +number of elements in the vectors. + +A further simplification of this formula, assuming the vectors X and Y being +row vectors, using matrix algebra (remember what you have learnt from matrix +multiplication): + + cov(X, Y) = X * Y' / (n-1) + +Sometimes we even forget about the normalization term `(n-1)`, so we end up +with `X*Y'`. + +Anyway, what we need to know about the covariance that it quantifies the extent +to which variables x and y tend to covary (what’s in a name?) linearly. In other +words, if the sign of the covariance is positive, it reflects the overall tendency +of y to go up if x goes up (and vice versa). If the sign of the covariance is +negative, y has the tendency to go down if x goes up, and vice versa. If the +covariance is very close to 0, x and y don’t tend to covary. + +The concept of covariance is very closely related to the concept of correlation +and regression, but this is not the time to discuss it in further detail. + +Let’s start with creating a few vectors in MATLAB to further illustrate the +concept of covariance. First, we will create two vectors x and y, with the mean +value subtracted: + + x = randn(1, 100); + x = x-mean(x); + + y = randn(1, 100); + y = y-mean(y); + +Now we can compute the covariance between x and y with just a simple command: + + covxy = (x*y')/99; + +We can also visualize these vectors: + + figure; hold on; plot(x); plot(y, 'r'); + figure; plot(x, y, '.'); + +The second figure is a so-called scatter plot, which shows that the individual +pairs of observations (the corresponding entries in the vectors x and y) are +distributed all over the Cartesian xy-plane. + +**_Q8.1.1 - Looking at the scatter plot with x along the horizontal and y along +the vertical axis. Are the two correlated?_** + +Now we can introduce covariance between the vectors x and y by adding some +shared signal to the vectors: + + a = 1; + b = 1; + c = randn(1, 100); + + x2 = x + a .* c; + x2 = x2 - mean(x2); + + y2 = y + b .* c; + y2 = y2 - mean(y2); + + covxy2 = (x2*y2')/99; + +**_Again make a scatter plot._** + +**_Q8.1.2 - What is the correlation coefficient between x2 and y2? Explain +the (approximate) number that you get (hint: think about how the correlation +relates to explained variance)._** + +Up until now we have explored the covariance between two variables represented +as vectors. It may come as no surprise that we can perform the same multiplication +operation to a pair of matrices. If we consider a sensor-level data matrix X +with the channels in the rows and the time points in the columns (and with the +mean across columns subtracted from each column), when we do `X*X'` we will +get an n x n covariance matrix (strictly speaking scaled with the number of +time points minus one: for the following let’s forget about this scaling, and +also let’s assume the individual rows to have a mean value of 0), where n is +the number of channels. Each element `c(i,j)` in this matrix, let’s call this +matrix `C` reflects the covariance between the i-th and j-th channel. + +There are a few remarks to be made about a covariance matrix: + +- When a covariance matrix is mirrored across the main diagonal of the matrix, you will get the same values. +- The elements on the main diagonal contain the auto-covariance between each channel and itself, and are always positive values. These values represent the variance of each channel. + +**_Q8.1.3 - Compute some normally distributed random data for 32 channels and +100 timepoints and represent this in a single matrix._** + +**_Q8.1.4 - Now compute the covariance between all channels; this can be represented +in a nchan x nchan matrix. You can use either for loops, or the corresponding +function in MATLAB._** + +**_Q8.1.5 - Plot the covariance matrix using `imagesc` and add a colorbar to +the figure. + +### 8.2 The sensor covariance gives information about the mixing process + +Let’s go back now to our basic observation model: `X = L*S+n`. If we now change +this equation to include the individual contributions of the single sources +s1, s2, ..., up to `sn`, and forget about the noise we get: + +X = l1*s1 + l2*s2 + ... + ln*sn + +When we compute `X*X'` and observe the individual terms in X we are essentially +doing the following: + +C = (l1*s1 + l2*s2 + ... +ln*sn) * (l1*s1 + l2*s2 + ... +ln*sn)' + +Note that here we are ignoring for now the normalization, i.e., the division +by N-1. + +There is a basic rule in matrix algebra that states that `(A*B)'` is the same +as `B'*A'`. Using this rule, and taking all the individual terms in the above +equation out of the brackets, we get: + + C = l1**s1**s1'*l1' + l2**s2**s1'*l1' + ... + ln**sn**s1'*l1' + ... + l1**s1**s2'*l2' + l2**s2**s2'*l2' + ... + ln**sn**s2'*l2' + ... + ... + l1**s1**sn'*ln' + l2**s2**sn'*ln' + ... + ln**sn**sn'*ln' + +Now the above equation looks complicated, but if we realize that `si*sj'` +represents the covariance (or auto-covariance) between the sources i and j, +and that these are scalar values (which means that it does not matter in which +order you multiply this in a matrix multiplication), we can rewrite the above +equation as: + + C = var(s1) **l1**l1' + cov(s2,s1)**l2**l1' + ... + cov(sn,s1)**ln**l1' + ... + cov(s1,s2)**l1**l2' + var(s2) **l2**l2' + ... + cov(sn,s2)**ln**l2' + ... + ... + cov(s1,sn)**l1**ln' + cov(s2,sn)**l2**ln' + ... + var(sn) **ln**ln' + +The terms `li*lj'` are also called vector outer-products, so in words the +above equation means that the covariance matrix can be represented as a _weighted +sum_ of the leadfield outer-products between all _pairs of sources_, where the +weights are based on the covariance between the corresponding pairs of sources. + +This is a very important feature of the covariance matrix that is exploited +by the beamformer, and is also exploited when we unmix the data matrix based +on a principal components analysis (which could be the topic of another session). + +### 8.3 Computing a beamformer solution + +After this section, you will + +- Understand the basic recipe for a beamformer analysis. +- Understand that the quality of the reconstruction critically depends on +the quality of the estimate of the covariance matrix. + +As already mentioned in the introduction section, the beamforming approach +is a so-called scanning method, where we investigate a set of locations in the +brain sequentially, where the question asked repeatedly is 'What would be the +best estimate of the dipole moment at _this_ location?’. + +For each of the locations investigated the beamformer algorithm computes a +spatial filter, also known as a set of weights (hence the convention to denote +the spatial filter with the letter w) that fulfills two criteria: + +1. The so-called filter gain is equal to 1. +2. The filter output should have minimum variance. + +The second criterion is specific to the type of beamformer we are typically +using, which is known as the minimum-variance beamformer. + +The first criterion of a filter gain of 1 is nothing else than 'what you get +is what you see’. In other words, the amplitude of the source signal that is +coming from a particular location in the brain (which is represented at the +channels by the leadfield from that location) should not be attenuated nor amplified +by the spatial filter that is estimated for that location. In a mathematical +formula this means `wi' * li = 1`. + +The second criterion aims that activity that is not coming from the location +under investigation does not affect ('leak into’) the estimate we’re currently +interested in. 'Filter output’ here means the estimated source activity. It +makes sense to require the _variance_ of this filter output to be as low as +possible, because if we are investigating a location from which no source activity +is emanating, we want the estimated activity to be as low as possible. + +It turns out that the two criteria expressed above are fulfilled by a spatial +filter computed as: + + w' = inv(lf'*inv(C)*lf) * lf' * inv(C) + +The beamformer makes an accurate reconstruction of the time course of the +underlying sources provided the individual sources are uncorrelated over time. +In reality neural sources are hardly even totally uncorrelated, but in practical +situations the beamformer is still tolerant against moderate levels of source +correlation (more about this later). It can be easily seen that the above formulation +satisfies criterion 1: + + w'*lf = inv(lf'*inv(C)*lf) * (lf'*inv(C)*lf) = I + +The proof that it also satisfies the second criterion can be demonstrated +with Lagrange multipliers (see https://en.wikipedia.org/wiki/Lagrange_multiplier), +but this falls beyond the scope of this course. + +The filter output or estimated source activity is defined as `w'*X` + +*Assume that we have data corresponding to 64 channels and and 1000 time points +and that we want to use the beamformer filter to estimate the activity of a +single dipole with known orientation, i.e., using a leadfield with 1 column.* + +**_Q8.3.1 - What is the size of the data matrix `X` and the size of the data +covariance matrix?_** + +**_Q8.3.2 - What is the size of the filter `w`?_** + +**_Q8.3.3 - What is the size of the result of the multiplication `w'*X`?_** + +### 8.4 Beamforming of simulated data. + +Let us now create some simulated data to demonstrate the beamformer. + + sensors = ni2_sensors('type', 'meg'); + headmodel = ni2_headmodel('type', 'spherical', 'nshell', 1); + + leadfield1 = ni2_leadfield(sensors, headmodel, [-2 -8 3 -1 0 0]); % position 1110 in the grid + leadfield2 = ni2_leadfield(sensors, headmodel, [-5 0 6 1 0 0]); % position 2342 in the grid + leadfield3 = ni2_leadfield(sensors, headmodel, [ 5 0 6 0 1 0]); % position 2352 in the grid + leadfield4 = ni2_leadfield(sensors, headmodel, [ 4 -2 7 0 0.2 0.7]); % position 2674 in the grid + [s1, t] = ni2_activation('latency', .45, 'frequency', 3); + [s2, t] = ni2_activation('latency', .50, 'frequency', 10); + [s3, t] = ni2_activation('latency', .50, 'frequency', 15); + [s4, t] = ni2_activation('latency', .55, 'frequency', 30); + + sensordata = leadfield1*s1 + leadfield2*s2 + leadfield3*s3 + leadfield4*s4 + randn(301, 1000)*1e-9; + +**_Q8.4.1 - Make a figure with 4 subplots underneath each other and plot the +activity of source 1 to 4 (along the vertical axis) versus time (along the horizontal +axis)._** + +**_Q8.4.2 - Make a figure with the sensordata versus time. Do you recognize +the source activity?_** + +We are first going to compute the beamformer spatial filters by hand, to demonstrate +the recipe. As mentioned in the previous section, the beamformer needs the sensor-level +covariance matrix and the leadfields for the locations at which the spatial +filters are going to be estimated. For now we will be using a 3-dimensional +regular 'grid’ of dipole positions, where the source locations are 1 cm apart. +Also, we will use FieldTrip to efficiently compute the leadfields. In principle +we could achieve the same using the `ni2_leadfield` function, but the first +approach is much faster: + + sourcemodel = ni2_sourcemodel('type', 'grid', 'resolution', 1); + + cfg = []; + cfg.sourcemodel = sourcemodel; + cfg.grad = sensors; + cfg.headmodel = headmodel; + sourcemodel = ft_prepare_leadfield(cfg); + L = cat(2, sourcemodel.leadfield{sourcemodel.inside}); + +**_Q8.4.3 - What is the size of the matrix L?_** + +We also need the covariance of the sensor-level data: + + M = mean(sensordata, 2); + M = repmat(M, 1, 1000); + C = ((sensordata-M)*(sensordata-M)')./999; + +Now we also compute the inverse of the covariance matrix, because it will +be used repeatedly. When it has to be recomputed each time again it slows down +the computations. + + iC = inv(C); + + for ii=1:size(L, 2)/3 + indx = (ii-1)*3 + (1:3); + Lr = L(:, indx); % Lr is the Nx3 leadfield matrix for source r + wbf(indx,:) = pinv(Lr'*iC*Lr)*Lr'*iC; + end + + sbf = wbf * sensordata; + +**_Q8.4.4 - Explain what the variable `indx` is and how it is used._** + +We can now inspect the reconstructed source time course at the locations at +which activity was simulated. Note that if we don't constrain the orientation +of the sources (i.e., use a 3-column leadfield per location) we will get a 3-row +spatial filter per location and hence also estimates for the x, y and z component +of the activity. + +Since we only did the computations for the grid positions inside the head, +we have to go from the original 3D grid indices that covers a square regular +grid (inside _and_ outside the head) to the list of estimations that are only +for the inside positions. + +**_Q8.4.5 - How many grid positions are there in the complete grid?_** + +**_Q8.4.6 - How many positions are there inside the head?_** + +**_Q8.4.7 - Make a plot of all grid positions and a plot of only the grid positions +that are inside the head. You can use plot3, axis equal and axis vis3d._** + +We now look up the grid indices at which the original sources were placed. + + index = 1:numel(sourcemodel.inside); % these are _all_ source positions, including the ones outside the brain + index1110 = find(index(sourcemodel.inside)==1110) % find the index of source position 1110, only considering the ones inside the brain + index2342 = find(index(sourcemodel.inside)==2342) + index2352 = find(index(sourcemodel.inside)==2352) + index2674 = find(index(sourcemodel.inside)==2674) + sel1110 = (index1110-1)*3 + (1:3) % find the three columns corresponding to source position 1110 + sel2342 = (index2342-1)*3 + (1:3) + sel2352 = (index2352-1)*3 + (1:3) + sel2674 = (index2674-1)*3 + (1:3) + +Each vector `selXXX` is a triplet of consecutive numbers that points to rows +in the matrix of wbf (and sbf) that we are going to explore first. + + figure + subplot(1, 2, 1); plot(t, sbf(sel1110,:)); + subplot(1, 2, 2); plot(t, s1); + figure + subplot(1, 2, 1); plot(t, sbf(sel2342,:)); + subplot(1, 2, 2); plot(t, s2); + figure + subplot(1, 2, 1); plot(t, sbf(sel2352,:)); + subplot(1, 2, 2); plot(t, s3); + figure + subplot(1, 2, 1); plot(t, sbf(sel2674,:)); + subplot(1, 2, 2); plot(t, s4); + +As you may have noticed, the resulting time courses are still a bit noisy +due to the noise in the data also being projected onto the estimated source +time courses. Since we only have 1000 samples, the structure of the noise in +the data is not so well estimated and hence cannot be so well suppressed. We +can improve the noise suppression using a regularized estimate of the inverse +covariance. We do that using Tikhonov regularization, also used in Ridge regression +(see https://en.wikipedia.org/wiki/Ridge_regression). For that we add a certain +amount of "identity matrix" to the covariance matrix prior to taking the inverse. +The amount of regularization can be determined from the plot of the eigenvalues +(or singular values) of the covariance, also known as the scree plot (see https://en.wikipedia.org/wiki/Scree_plot). + + iCr = inv(C + 1e-18*eye(301)); + +### 8.5 Beamformer and depth bias + +After this section, you will + +- Understand the depth bias of the beamformer, and its cause. +- Understand how in the typical application a contrast between two conditions is needed for a meaningful source reconstruction. +- Understand how the performance of the beamformer is affected by correlation of the source time courses. + +In the previous section we visualized the reconstructed time courses at the +locations close to where we actually simulated the activity. In real applications, +these locations of course are unknown, and we typically need to look for local +maxima in the reconstructed 'image’. The quantity that is used is the filter +output, i.e., the variance of the source time courses, which is also often referred +to as the 'power’ of the sources. Apart from the fact that in real applications +the data are usually quite noisy, causing difficulties in the estimation of +the local maxima, there is an important feature of the beamformer, which further +complicates the evaluation of peaks in the image. This is the so-called depth +bias. + +To illustrate this depth bias, we first restructure the source-reconstructed +simulated data of section 3.2 in a data-structure that can be used by FieldTrip, +so that we can use FieldTrip's visualization functionality. If you don’t have +the `sbf` and `sourcemodel` variables in MATLAB memory anymore, you need to +reconstruct these in order to be able to proceed. + +In order to represent the reconstruction as an image, we need to express the +variance of the sources at each of the locations in a single number. Recall +that at each location of the 3-dimensional grid the source activity consists +of three time courses, one for each 'cardinal’ x/y/z orientation. A common way +to achieve this is to sum the variance across the three orientations. This is +essentially the application of Pythagoras’ rule (without explicitly taking the +square and the square root, since variance is already a 'squared’ value). With +our sbf variable it is possible to do this in the following way. Note that this +variable is a 'number of inside sources x 3’ times 'number of timepoints’ matrix. +If we compute the variance across time (var(sbf, [], 2)) we end up with a vector, +that in consecutive triplets contains the variance in the x/y/z orientation +at the 'inside’ dipole locations of the sourcemodel. We can now efficiently +create the variance for each location by using MATLAB’s reshape and sum functions: + + pbf = var(sbf, [], 2); + pbf = reshape(pbf, 3, numel(pbf)/3); + pbf = sum(pbf, 1); + +Take a moment to try and understand what is going on in the second and third +line. + +Now, create a FieldTrip structure of the 'source’ type, and use `ft_sourceplot` +to visualize this. + + source = []; + source.pos = sourcemodel.pos; + source.dim = sourcemodel.dim; + source.inside = sourcemodel.inside; + source.inside = reshape(source.inside, source.dim); + source.pow = zeros(size(source.pos, 1), 1); + source.pow(source.inside) = pbf; + source.pow = reshape(source.pow, source.dim); + + cfg = []; + cfg.funparameter = 'pow'; + cfg.method = 'slice'; + cfg.nslices = 10; + cfg.funcolorlim = [0 0.2]; + ft_sourceplot(cfg, source); + +**_Q8.5.1 - Include the figure that you just created in your answers._** + +What you see now is a set of slices through the spherical volume conductor, +where each slice is parallel to the xy-plane. It starts at the left-lower bottom +with the slice that has the lowest z-coordinate, and ends at the upper right +with the crown of our 'simulated’ head. The x-axis is running from left to right, +and the y-axis is running from bottom to top. + +The most prominent feature in this image is the enormous 'blob’ in the middle +of the 'head’. This represents the depth bias of the beamformer, and is an important +motivation for why most successful applications of the beamformer always involve +the reconstruction of a contrast between two conditions (another important motivation +is that we are usually interested in the features in our data that change as +a consequence of an experimental manipulation, which by construction always +involves the comparison across two (or more) experimental conditions). In this +simulated data example with a relatively high signal-to-noise ratio of the data, +admittedly three out of the four simulated sources are visible. + +The depth bias is directly caused by the unit gain constraint of the spatial +filter. The further away the source is from the sensors, the smaller the numbers +in the leadfield (why?). If the dot-product between the spatial filter and the +leadfield is required to be 1, the small values in a deep source’s leadfield +need to be counterbalanced by large values in the spatial filter. As a direct +consequence the spatial filter’s output (`w'**C**w`) will be much larger for a +deep source than for a more superficial source. + +### 8.6 Accounting for depth bias: computing a contrast + +In the previous section we have seen that the depth bias is related to the +fact that the spatial filters for deep sources have much higher values than +the spatial filters for superficial sources. The depth bias is therefore a feature +of the spatial filters. One way to account for this is to normalize the estimated +source power by an estimate of how noise projects through the spatial filters. +This normalized quantity is known as the 'neural activity index’, but in real +data applications this usually does not work that well. The reason is that a +good estimate of the noise is difficult to obtain, and usually the noise covariance +is assumed to be a scaled identity matrix (this reflects noise to be uncorrelated +across channels, with the same variance for each of the channels), which is +far from reality. A better strategy is to compute a contrast between two conditions. +The idea is that the features in the data that are the same across conditions +(such as the depth bias) will cancel each other when making the contrast. + +We will use the variables `sensordata`, `sourcemodel` and `L` that we also +used in section 3. We also need the leadfields and source time courses that we +used for the simulations. If you don’t have these variables anymore in MATLAB +memory, you should re-create them. + +Now we will simulate data from a 'second’ condition (compared to the original +variable sensordata), where the sources have the exact same locations and time +courses, but the amplitude of two sources is decreased, and the amplitude of +the two other sources is increased, relative to the 'first’ condition. + + sensordata2 = 1.25 .* leadfield1*s1 + ... + 0.80 .* leadfield2*s2 + ... + 0.80 .* leadfield3*s3 + ... + 1.25 .* leadfield4*s4 + ... + randn(301, 1000)*1e-9; + +We will now compute the spatial filters using the covariance estimated from +the data combined across the two conditions. In this way we will end up with +a single set of spatial filters, and thus, in comparing across the two conditions, +we ensure that the depth bias is the same for condition 1 and 2. The covariance +of the data combined across the 2 conditions can be computed in several ways, +e.g. by averaging the single condition covariances. Here, we first concatenate +the data and subsequently compute the covariance. Now we will use the MATLAB +cov function, rather than computing the covariance 'by hand’ (i.e., by first +computing the mean across time etc.). + + C2 = cov([sensordata sensordata2]'); + iCr2 = inv(C2 + eye(301)*2e-18); + + for ii=1:size(L, 2)/3 + indx = (ii-1)*3 + (1:3); + Lr = L(:, indx); % Lr is the Nx3 leadfield matrix for source r + wbf(indx,:) = pinv(Lr'*iCr2*Lr)*Lr'*iCr2; + end + + sbf1 = wbf * sensordata; + sbf2 = wbf * sensordata2; + + pbf1 = var(sbf1, [], 2); + pbf1 = reshape(pbf1, 3, numel(pbf1)/3); + pbf1 = sum(pbf1, 1); + + pbf2 = var(sbf2, [], 2); + pbf2 = reshape(pbf2, 3, numel(pbf2)/3); + pbf2 = sum(pbf2, 1); + +We can now create a FieldTrip style source-structure, and store in the field +pow a measure of the difference between condition 1 and 2. + + source = []; + source.pos = sourcemodel.pos; + source.dim = sourcemodel.dim; + source.inside = sourcemodel.inside; + source.inside = reshape(source.inside, source.dim); + source.pow = zeros(size(source.pos, 1), 1); + source.pow(source.inside) = (pbf1-pbf2)./((pbf1+pbf2)/2); + source.pow = reshape(source.pow, source.dim); + + cfg = []; + cfg.funparameter = 'pow'; + cfg.method = 'slice'; + cfg.nslices = 10; + cfg.funcolorlim = [-0.3 0.3]; + ft_sourceplot(cfg, source); + +**_Q8.6.1 - Include the figure that you just created in your answers._** + +### 8.7 Beamforming and correlated sources + +An unfortunate characteristic feature of the beamformer is its inability to +deal with strongly correlated sources. This is a direct consequence of the fact +that, in the presence of correlated sources, the leadfield outer products of +the sources that are correlated contribute significantly to the sensor level +covariance matrix (see section 2.2). As a result, the spatial filters, which +aim to suppress interfering (but uncorrelated) activity originating from other +locations lead to a distortion of the true source time courses. In this section +we will start from a fresh simulation, so you’d best issue a `clear all` command +before you proceed. + +First, we create some sensor data: + + sensors = ni2_sensors('type', 'meg'); + headmodel = ni2_headmodel('type', 'spherical', 'nshell', 1); + leadfield1 = ni2_leadfield(sensors, headmodel, [-5.3 0 5.9 1 0 0]); % position 2342 in grid + leadfield2 = ni2_leadfield(sensors, headmodel, [4.9 0 6.2 0 1 0]); % position 2352 in grid + +create the time course of activation + [s1, t] = ni2_activation('latency', .5, 'frequency', 10); + [s2, t] = ni2_activation('latency', .49, 'frequency', 10); + +create the sensor data + sensordata = leadfield1*s1 + ... + leadfield2*s2 + ... + randn(301, 1000)*0.04e-9; + +Now we can assess the correlation between s1 and s2, and we will note that +it is significant: + + c = corr(s1', s2'); + +**_Q8.7.1 - What is the correlation?_** + +**_Q8.7.2 - Plot the source time course s1 and s2 in a single figure._** + +We can make a scree plot of the eigenvalues of the sensor covariance. It reveals +that there is largely only one source component contributing to the data + + v = svd(cov(sensordata')); + figure; plot(v(1:10), 'o') + +**_Q8.7.3 - Try making sensordata with more and with less correlated sources._**