OPM tutorial for CuttingEEG X

workshop/cuttingeegx.md

@@ -9,7 +9,142 @@ CuttingEEG is turning 10 years old! This milestone calls for a special edition,
 
 Alongside the [CuttingEEG X](https://cuttingeegx.org) conference we will be organizing some local workshops in the morning.
 
-- Who: Julia Chauvet, Konstantinos Tsilimparis, Britta Westner, Tilman Stephani, Robert Oostenveld
-- When: 30 to 31 October 2024
+- Who: Konstantinos Tsilimparis, Robert Oostenveld
+- When: 31 October 2024
 - Where: Nijmegen
 - See <https://cuttingeegx.org/registration/#Nijmegen> for more details
+
+## The tutorial
+
+We will run a tutorial on [preprocessing and coregistration of SQUID-based and OPM-based data](/workshop/cuttingeegx/squids_vs_opms). We will explore the differences in data analysis between the two systems. 
+
+## Getting started with the hands-on session
+
+In this workshop you will work on your own laptop computer. It will be an in-person event with no possibilities for hybrid or online attendance.
+
+### Wifi access
+
+If you need wifi access and you don't have a eduroam account through your institution, it is possible to get a visitor access. This needs to be renewed each day. Please follow the instructions on [this intranet page](https://intranet.donders.ru.nl/index.php?id=eva).
+
+### MATLAB
+
+For the hands-on sessions we assume that you have a computer with a relatively recent version of MATLAB installed (preferably < 5 years old, >= 2019a/b). 
+
+### FieldTrip
+
+To get the most recent copy of FieldTrip, you can follow this [link](https://github.com/fieldtrip/fieldtrip/releases/tag/20240916), download the zip-file, and unzip it at a convenient location on your laptop's hard drive. 
+
+Alternatively, you could do the following in the MATLAB command window. 
+
+```
+% create a folder that will contain the code and the data, and change directory
+mkdir('cuttingeegx');
+cd('cuttingeegx');
+
+% download and unzip fieldtrip into the newly created folder
+url_fieldtrip = 'https://github.com/fieldtrip/fieldtrip/releases/tag/20240916.zip';
+unzip(url_fieldtrip);
+```
+
+Upon completion of this step, the folder structure should look something like this: 
+
+```bash
+cuttingeegx/
+|-- fieldtrip-20240916
+    |-- bin
+    |-- compat
+    |-- connectivity
+    |-- contrib
+    |-- external
+    |-- fileio
+    |-- forward
+    |-- inverse
+    |-- plotting
+    |-- preproc
+    |-- private
+    |-- qsub
+    |-- realtime
+    |-- specest
+    |-- src
+    |-- statfun
+    |-- template
+    |-- test
+    |-- trialfun
+    |-- utilities
+```
+
+{% include markup/red %}
+If you have downloaded and unzipped by hand, it could be that there's an 'extra folder layer' in your directory structure. We recommend that you remove this extra layer, i.e. move all content one level up.
+{% include markup/end %}
+
+### Test your installation in advance
+
+To have a smooth experience - and to avoid having to spend precious debugging time during the hands-on sessions - we recommend that you [test your MATLAB and FieldTrip installation in advance](/workshop/toolkit2024/test_installation).
+
+## The data used in this tutorial
+
+Next, we proceed with downloading the relevant data. The data that are used in the hands-on sessions, are stored on the FieldTrip [download-server](https://download.fieldtriptoolbox.org/workshop/cuttingeegx). You can either ‘click around’ using web browsers and/or explorer windows to grab the data that are needed, or instead (less work, at least if it works) execute the MATLAB code below.
+
+Please ensure that your present working directory is the ```cuttingeegx``` folder, which you created in the previous step.
+
+```
+% create a folder (within cuttingeegx) that will contain the data, to keep a clean structure
+mkdir('data');
+cd('data');
+
+% create a folder and download the SQUID and OPM dataset
+url_tutorial = 'https://download.fieldtriptoolbox.org/workshop/cuttingeegx';
+fnames = {''};
+for k = 1:numel(fnames)
+  websave(fnames{k}, fullfile(url_tutorial, fnames{k}));
+end
+cd('../');
+```
+
+At this stage, you ideally have a directory structure that looks like the following one:
+
+```bash
+cuttingeegx/
+|-- data
+    |-- opm
+    |-- squid
+    |-- mri
+|-- fieldtrip-20240916
+    |-- bin
+    |-- compat
+    |-- connectivity
+    |-- contrib
+    |-- external
+    |-- fileio
+    |-- forward
+    |-- inverse
+    |-- plotting
+    |-- preproc
+    |-- private
+    |-- qsub
+    |-- realtime
+    |-- specest
+    |-- src
+    |-- statfun
+    |-- template
+    |-- test
+    |-- trialfun
+    |-- utilities
+```
+
+
+Whenever starting a fresh MATLAB session, to configure the right FieldTrip paths, execute the following: 
+
+```
+% change into the 'cuttingeegx' folder and then do the following
+restoredefaultpath
+addpath('fieldtrip-20240417');
+addpath(genpath('data'));
+ft_defaults;
+```
+
+ `ft_defaults` command ensures that all of FieldTrip's required subdirectories are added to the path.
+
+{% include markup/red %}
+Furthermore, please do NOT add FieldTrip with all subdirectories, subdirectories will be added automatically when needed, and only when needed (see this [FAQ](/faq/installation).
+{% include markup/end %}

workshop/cuttingeegx/squids_vs_opms.md

@@ -5,4 +5,238 @@ tags: [cuttingeegx]
 
 # SQUIDs versus OPMs
 
-To be implemented, see https://github.com/fieldtrip/website/issues/798
+## Introduction
+
+Traditional MEG systems use superconducting quantum interference devices (SQUIDs), which require expensive cryogenic cooling with liquid helium and involve placing the sensors 1.5-2 cm from the scalp. Optically pumped magnetometers (OPMs), a new generation of MEG, eliminate the need for cryogenic cooling, allowing sensors to be placed closer to the scalp. OPMs are small individual sensors that allow the researcher to choose where to place them over the head. These sensors are held in place by 3D-printed helmets, which can be customized for each subject.
+
+Even though, both SQUID and OPM system detect the brain's magnetic fields, their data analysis steps differ a bit. In this tutorial we explore the differences in preprocessing and coregistration between the two systems.
+
+OPMs are flexible in their placement allowing for new recording strategies. These new recording strategies lead to new data analysis strategies. For example, in this tutorial we use a small number of OPMs (32 sensors) and do sequential recordings in which we position them over different places over the scalp.
+
+OPMs are magnetometers, as their name suggests. Magnetometers are more sensitive to environmental noise than gradiometers, which most SQUID systems have. Several data analysis algorithms to remove environmental noise have been proposed (see [Seymour et al (2022)](https://www.sciencedirect.com/science/article/pii/S1053811921011058?via%3Dihub) for more details). In this tutorial, we apply homogeneous field correction (HFC). HFC works better with a large number of sensors (more than 32).
+
+Unlike SQUID systems, which have standard coregistration procedures, OPMs don't have a single standard. In this tutorial, we coregister the OPMs with the MRI using an optical 3D scanner which captures the participant’s facial features along with the OPM helmet (Zetter et al., 2019).
+
+
+This tutorial does not cover follow-up analyses (like source reconstruction) which in principle should not differ from the SQUID follow-up analyses, or alternative coregistration methods which are covered in the tutorial on [coregistration of Optically Pumped Magnetometer (OPM) data](tutorial/coregistration_opm/).
+
+## Background
+
+
+In this tutorial we will use recordings made with 32 OPM sensors  placed in an adult-sized “smart” helmet with a total of 144 slots. Each slot allows the sensor to move along a single axial direction. That way, the sensor can slide in its slot until it touches the head surface, regardless of the head size and shape. To limit head movements we mounted the helmet on a wooden plate.
+
+To acquire a measurement for each of the 144 helmet slots, we divided the experiment into six sequential recordings while maintaining the participant's head in a fixed position. We kept 9 sensors around the participant’s head fixed for each recordings. Since the OPM sensors touched the participant’s head and the helmet was mounted, these 9 sensors were able to keep the participant's head fixed throughout the experiment. In each recording we moved the remaining 23 sensors to fill every helmet slot.
+
+### The dataset used in this tutorial
+The data for this tutorial was recorded with a 32-sensor FieldLine HEDscan v3 system with a so-called smart helmet. Each OPM sensor has one channel that measures the normal component of the magnetic field. 
+
+We perform an experiment with left median nerve stimulation on a single participant using both the SQUID and the OPM system. We expect the activity to be modelled 20 ms post-stimulation with a dipole located at the right Primary Somatosensory (S1) area (Andersen & Dalal, 2021; Boto et al., 2017; Buchner et al., 1994)). 
+
+The dataset can be downloaded from our download server.
+
+## Procedure 
+In this tutorial for the SQUIDs we will take the following steps:
+- Define trials and read the data using **[ft_definetrial](/reference/ft_definetrial)** and **[ft_preprocessing](/reference/ft_preprocessing)**
+- Removing artifacts using **[ft_rejectvisual](/reference/ft_rejectvisual)**
+- Compute the averaged ERFs using **[ft_timelockanalysis](/reference/ft_timelockanalysis)**
+- Visualize the results for all the channels with **[ft_multiplotER](/reference/ft_multiplotER)**
+- Plot the 2D sensor topography for a specified latency with **[ft_topoplotER](/reference/ft_topoplotER)**
+- Coregister MRI with SQUIDs using **[ft_volumerealign](/reference/ft_volumerealign)**
+- Plot the 3D sensor topography for a specified latency with **[ft_plot_topo3d](/reference/plotting/ft_plot_topo3d)**
+
+For the OPMs we will take the following steps:
+- Define trials and read the data using **[ft_definetrial](/reference/ft_definetrial)** and **[ft_preprocessing](/reference/ft_preprocessing)**
+- Removing artifacts using **[ft_denoise_hfc](/reference/ft_denoise_hfc)** and **[ft_rejectvisual](/reference/ft_rejectvisual)**
+- Compute the averaged ERFs using **[ft_timelockanalysis](/reference/ft_timelockanalysis)**
+- Renaming duplicate channels
+- Concatenate the data over the six recordings using **[ft_appendtimelock](/reference/ft_appendtimelock)**
+- Add NaNs to the missing channels
+- Visualize the results for all the channels with **[ft_multiplotER](/reference/ft_multiplotER)**
+- Plot the 2D sensor topography for a specified latency with **[ft_topoplotER](/reference/ft_topoplotER)**
+- Coregister MRI with OPMs using an optical 3D scanner. For this several functions are used: **[ft_volumerealign](/reference/ft_volumerealign)**, **[ft_read_headshape](/reference/fileio/ft_read_headshape)**, **[ft_meshrealign](/reference/ft_meshrealign)**, **[ft_defacemesh](/reference/ft_defacemesh)**, and **[ft_transform_geometry](/reference/utilities/ft_transform_geometry)**.
+- Append sensors from the six recordings using **[ft_appendsens](/reference/ft_appendsens)**
+- Plot the 3D sensor topography for a specified latency with **[ft_plot_topo3d](/reference/plotting/ft_plot_topo3d)**
+
+
+## SQUID
+### Preprocessing & trial definition
+
+We begin by loading the SQUID data and defining trials. In the experiment, the inter-trial interval ranged from 800-1200 ms. We select a 200 ms prestimulus and 400 ms poststimulus window. We then select trials where left median nerve stimulation occurred (trigger code = 1).
+
+```
+cfg                     = [];
+cfg.dataset             = 'M:\Documents\cuttingeegx-workshop\data\squid\subj001ses002_3031000.01_20231208_01.ds';
+cfg.trialdef.eventtype  = 'UPPT001';
+cfg.trialdef.eventvalue = [1]; % select stimulation trials only (which are ~85 % of all the trials)
+cfg.trialdef.prestim    = 0.2; % ISI varied from 0.8 to 1.2 sec
+cfg.trialdef.poststim   = 0.4;
+cfg                     = ft_definetrial(cfg);
+
+cfg.demean         = 'yes';  
+cfg.baselinewindow = [-Inf 0]; 
+cfg.continuous     = 'yes'; % without that I get the error: 'requested data segment extends over a discontinuous trial boundary'. I guess I need cfg.continuous since 'subj001ses002_3031000.01_20231208_01.ds' is initially a continuous recording and I do not want to lose the ITI during preprocessing.
+data_stim          = ft_preprocessing(cfg);
+
+save data_stim data_stim
+```
+
+### Removing artifacts
+
+The summary method in **[ft_rejectvisual](/reference/ft_rejectvisual)** allows to manually remove trials and/or channels that have high variance. 
+
+```
+load data_stim data_stim
+
+cfg             = [];
+cfg.method      = 'summary';
+data_stim_clean = ft_rejectvisual(cfg,data_stim);
+
+save data_stim_clean data_stim_clean
+```
+
+### Compute ERFs
+
+We start by removing trials that have higher variance than 8e-25 Tesla-squared, then assess the channels. Removing trials also reduces channel variance, so no need to reject any channels. We then calculate the ERFs: 
+
+```
+load data_stim_clean data_stim_clean
+
+cfg     = [];
+avg_stim = ft_timelockanalysis(cfg, data_stim_clean);
+
+save avg_stim avg_stim
+```
+### Visualisation
+
+We plot the activity from all the sensors.
+
+```
+load avg_stim avg_stim
+
+% find the time window of the activity
+cfg = [];
+cfg.layout = 'CTF275_helmet';
+ft_multiplotER(cfg, avg_stim); % use interactive
+```
+
+We are now going to use the interactive feature of **[ft_multiplotER](/reference/ft_multiplotER)** to find our activity of interest. We select the sensors that are on top of the right primary somatosensory area since that is where we expect our activity to be localised. We see a negative peak around 20 ms (more specifically 35-50 ms) post-stimulation. We can select this time window to see the topography. The dipolar pattern the right primary somatosensory area is now visible. We can also plot this dipolar pattern with **[ft_topoplotER](/reference/ft_topoplotER)**.
+
+```
+% plot the activity at [0.035 0.050]
+cfg = [];
+cfg.xlim = [0.035 0.050];
+cfg.layout = 'CTF275_helmet';
+ft_topoplotER(cfg, avg_stim); 
+
+print -dpng figures/cuttingeegx_topo_squid.png
+```
+
+### Coregister the anatomical MRI to the SQUID coordinate system
+
+We read the MRI, SQUID channels and the Polhemus headshape in memory. I recommend converting all quantities to SI units to ensure consistent units throughout your pipeline.
+
+```
+mri    = ft_read_mri('M:\Documents\cuttingeegx-workshop\data\dicoms\00001_1.3.12.2.1107.5.2.19.45416.2022110716263882460203497.IMA');
+ctf275 = ft_read_sens('M:\Documents\cuttingeegx-workshop\data\squid\subj001ses002_3031000.01_20231208_01.ds', 'senstype', 'meg');
+shape  = ft_read_headshape('M:\Documents\cuttingeegx-workshop\data\squid\subj001ses001ses002.pos');
+
+mri = ft_convert_units(mri, 'm');
+ctf275 = ft_convert_units(ctf275, 'm');
+shape  = ft_convert_units(shape, 'm');
+
+save mri mri
+save ctf275 ctf275
+save shape shape
+```
+
+We coregister the MRI with the SQUIDs by converting the MRI coordinate system to match that of the SQUIDs. In other words, we ensure that the three fiducials (nasion, LPA and RPA) defining the MRI coordinate system are aligned to the same three fiducials defining the SQUID coordinate system.
+
+```
+load shape shape
+load ctf275 ctf275
+
+% assign the 'ctf' coordinate system to the MRI 
+cfg          = [];
+cfg.method   = 'interactive';
+cfg.coordsys = 'ctf';
+mri_realigned1 = ft_volumerealign(cfg, mri);
+
+save mri_realigned1 mri_realigned1
+```
+
+
+To refine the coregistration, we can also coregister the Polhemus head shape with the skin surface that is extracted from the previously aligned MRI. Now not only the fiducials but also other skin surface points acquired with the Polhemus are used for the alignment.
+
+```
+load mri_realigned1 mri_realigned1
+
+cfg = [];
+cfg.method = 'headshape';
+cfg.headshape.headshape = shape;
+mri_realigned2 = ft_volumerealign(cfg, mri_realigned1); 
+
+% check 
+cfg = [];
+cfg.method = 'headshape';
+cfg.headshape.headshape = shape;
+cfg.headshape.icp = 'no';
+ft_volumerealign(cfg, mri_realigned2); % takes too long
+
+save mri_realigned2 mri_realigned2
+```
+
+
+### 3D sensor topography
+Let's plot the sensor topography in 3D using **[ft_plot_topo3d](/reference/plotting/ft_plot_topo3d)**. This will help us visualize where the sensors are positioned relative to the head. To make this clearer, we'll also display the scalp and brain surfaces.
+
+First, we'll prepare the surface mesh of the scalp and brain:
+
+```
+load mri_realigned2 mri_realigned2
+
+cfg                = [];
+cfg.output         = {'brain', 'scalp'};
+mri_segmented          = ft_volumesegment(cfg, mri_realigned2);
+
+cfg                = [];
+cfg.tissue         = {'brain'};
+mesh_brain         = ft_prepare_mesh(cfg, mri_segmented);
+
+cfg                = [];
+cfg.tissue         = {'scalp'};
+mesh_scalp         = ft_prepare_mesh(cfg, mri_segmented);
+
+save mesh_brain mesh_brain
+save mesh_scalp mesh_scalp
+```
+Now, we'll plot the 3D sensor topography for the time window [0.035, 0.050] seconds, along with the brain and scalp:
+
+```
+sampling_rate = 1200; % in Hz
+prestim = 0.2;
+
+I1 = (prestim+0.035)*sampling_rate; 
+I2 = (prestim+0.050)*sampling_rate;
+
+selected_avg = mean(avg_stim.avg(:, I1:I2), 2);
+
+figure;
+ft_plot_mesh(mesh_scalp, 'facealpha', 0.5, 'facecolor', 'skin', 'edgecolor', 'none', 'edgecolor', 'skin' )
+hold on
+ft_plot_mesh(mesh_brain, 'facecolor', 'brain', 'edgecolor', 'none');
+hold on
+ft_plot_topo3d(pos275,selected_avg, 'facealpha', 0.9)
+camlight
+view([360 0])
+
+
+print -dpng M:\Documents\cuttingeegx-workshop\code\figures\squid\cuttingeegx_topo3d_squid.png
+```
+
+The SQUID sensors are located 1.5-2.0 cm from the scalp, resulting in a less focal sensor topography.
+
+## OPM
+### Preprocessing
+HFC is a method for denoising MEG data based on a spatially homogeneous model of the background magnetic field across the OPM array. This method has previously been used successfully for reducing magnetic interference in OPM magnetometers (see Supplementary Material Fig. A3; Hill et al., 2022; Mellor et al., 2023; Seymour et al., 2022) that are more sensitive to the environmental noise than the SQUID gradiometers. 
+
+### Coregistration