Pipelines Overview

Different pipelines can be built using the available software at NYUAD MEG Lab tailored to best answer your research questions.

Important

Pipelines provided for use knowing that bugs could exist

Provided pipelines by the MEG lab can never be tested for all situations and cases, therefore it is for the project owner’s responsability to check whether or not the output of any used pipeline is correct. The pipelines users should signal any suspiscious output, or suggestion about test cases to make the test benchmarks more robust by raising an issue on the repository.

Software stack and installation

The folowing software/library are available MEG/EEG data analysis:

  • MNE Python library

  • FieldTrip

  • BESA

Sample pipelines are provided for each one of them

MNE Python

Installation

The MNE-based pipelines in this repository run in a dedicated Conda environment named neurowaves, defined by pipeline/mne_pipelines/kit_general_pipelines/environment.yml. Creating, verifying and updating it is described on KIT general pipeline: environment setup and running. The MNE standalone installer is no longer the recommended way to run the lab pipelines: it lacks some of the pipeline’s dependencies and its versions cannot be pinned by the repository.

FreeSurfer for source localisation

Source localisation with MNE needs a FreeSurfer reconstruction of each participant’s T1w MRI. FreeSurfer runs on Linux (on Windows, inside WSL2), and stage 7 of the KIT general pipeline runs it for a whole BIDS dataset, and stage 7-a has every reconstruction checked by a reviewer before the MEG-MRI coregistration (stage 8) uses it. Its installation is described on Optional: FreeSurfer for stage 7.

Starting a project with the general MNE-KIT pipeline

You have acquired MEG-KIT data at NYUAD, organised it in BIDS following the lab’s naming convention (see Data Naming and Storing: directory structure, entity order, and the naming rules for .con, .mrk, laser-scan head shape and empty-room files), and want to start pre-processing. The pipeline finds every file through these BIDS names, so a dataset that does not follow the convention will be partly or entirely skipped. The layout it expects is:

<datasets parent>/                 # value of MEG_DATA or project.root_override
    <project name>/                # project.name: the BIDS root
        dataset_description.json
        participants.tsv
        sub-001/
            anat/sub-001_T1w.nii.gz        # only needed for stage 7
            ses-01/                        # only if the subject has several sessions
                meg/
                    sub-001_ses-01_task-<task>_run-01_meg.con
                    sub-001_ses-01_task-<task>_run-01_proc-CALMnoisereduction_meg.con
                    sub-001_ses-01_task-rest_acq-<N>_space-ALS_markers.mrk
                    sub-001_ses-01_acq-head_headshape.txt
                    sub-001_ses-01_acq-points_headshape.txt
        derivatives/               # created by the pipeline, one folder per stage

The KIT general pipeline is a sequence of numbered scripts driven by one YAML config per dataset, passed to every stage with -c / --config. Before running any stage, create a preliminary config for your project as follows.

  1. Copy the template and name it after your project:

    cd pipeline/mne_pipelines/kit_general_pipelines/pipeline_config_files
cp config_template.yml <project>_config.yml

    Relative paths written inside the YAML are resolved from the folder that contains the YAML, not from the folder you launch the scripts from.

  2. Point it at your dataset (project block). The pipeline builds the BIDS root as <parent>/<project.name>:

    • name: the exact name of your BIDS root folder (case-sensitive on Linux), for example "kidlang".

    • root_env / root_override: where <parent> comes from. Either keep root_override: null and set the environment variable named by root_env (default MEG_DATA) to the parent folder in the shell that runs the stages (or before starting VS Code):

      export MEG_DATA=/path/to/datasets_parent      # Linux / macOS
$env:MEG_DATA = "D:\MEG_Data"                  # Windows PowerShell

      or write the parent folder directly in root_override. On Windows use forward slashes, single quotes or doubled backslashes ("D:/MEG_Data", 'D:\MEG_Data' or "D:\\MEG_Data"): a single backslash inside double quotes breaks the YAML.

  3. Choose the subjects (subjects block). Labels are written without the sub- prefix. include: [] processes every sub-* folder found in the BIDS root; exclude always skips the listed labels (pilots, unusable subjects). For a first run, restrict include to one or two subjects, for example include: ["001"].

  4. Choose the recordings (bids_selection block). Every filter left as [] means “all”. Values are the labels as they appear in the file names, without the entity prefix ("01", not "ses-01"). Set at least:

    • tasks: your task label(s) if the dataset holds more than one task.

    • processings: whether to process the raw or the CALM noise-reduced .con files. [null] selects raw files only (no proc- entity), [CALMnoisereduction] the CALM files only, and [] processes both as separate recordings.

    The other filters (sessions, acquisitions, runs, splits, recordings, spaces, descriptions) can stay empty unless you want to process a subset. The nested exclude block removes matches, for example exclude: {tasks: [rest]}.

  5. Reports (reports block). Set formats: ["html"] unless PDF reports are set up (see Optional: PDF reports).

  6. Leave every other block at its template value for now (kit2fiff, bad_channels, filters, events, epoching, plots, freesurfer). Do not delete them: they are read by later stages.

The keys changed in a typical preliminary config are therefore only:

project:
  name: "myproject"
  root_env: "MEG_DATA"
  root_override: null
subjects:
  include: ["001"]
  exclude: []
bids_selection:
  tasks: ["mytask"]
  processings: [null]
reports:
  formats: ["html"]

Every stage prints the BIDS root it resolved when it starts. On your first run, check that this path is your dataset folder; if it is not, fix project.name or the parent folder before going further.

The settings that are specific to each processing step are described on the page of that step below, and you complete the same YAML as you progress:

  • kit2fiff (stage 2): trigger decoding (stim, slope, stimthresh, stim_code), coreg.mode: automatic | manual and n_jobs for the automatic mode.

  • events.multi_channel.max_onset_skew_ms (stages 1 and 5): grouping of binary-coded triggers.

  • bad_channels (stage 3-a): segment plan, detector thresholds (flat/peak, variance, PSD, line noise, LOF, Maxwell), decision rule (decision.min_detectors), per-subject manual and ignore_detected lists, and the review browser settings.

  • filters (stage 4): notch and band-pass. When the block is absent, stage 4 applies a notch at the 50 Hz harmonics and a 1 to 40 Hz band-pass.

  • events and plots (stage 5): the desc label written by stage 5 and read by stage 6, the source branch used by stage 6, find_events kwargs for the STI 014 variant, and the stim-channel plots.

  • epoching (stage 6): global tmin, tmax and baseline plus per-event-code windows.

  • freesurfer (stage 7): missing_anat: fsaverage | error, resource planning (jobs, threads_per_subject), auto_cw256 and the mkheadsurf scalp surface (scalp_surface, scalp_threshold).

  • freesurfer_qc (stage 7-a): flag thresholds, figures, freeview button and YAML verdict overrides.

  • coregistration (stage 8): require_mri_qc and the QA figures.

  • reports (all stages): per-script switches, quarto_env and output_dir.

The commented config_template.yml remains the reference for every key, and kidlang_config.yml is a complete worked example.

Order of pages for the general MNE-KIT pipeline

Read the pages in this order:

  • Environment setup and running: create the neurowaves environment, run the stages, find the outputs and the auto-generated report.

  • Stage 2, KIT to FIFF: convert the raw .con files (with and without CALM noise reduction) to .fif with head shape and markers, automatically or through the GUI.

  • Stage 1, trigger count sanity check (single-channel and binary-coded trigger modes): verify that no triggers were lost and that they are detected in the right order. Recommended before any analysis of analog-trigger datasets.

  • Stage 3-a, bad channel detection: human-in-the-loop audit of the 207 KIT gradiometers (inspection in MNE’s browser, automatic candidates, confirmation); the resulting list is applied by stage 4 and inherited by every later stage.

  • Stage 3-b, ICA: COGNESTIC-style ICA (fit on a 1 Hz high-passed copy, ECG/EOG proposals, review in MNE’s source browser); stage 4 applies the reviewed exclusions to the unfiltered data before filtering.

  • Stage 5, compute events: write the .eve files from the analog trigger lines or from the STI 014 channel; stage 6 epochs from them.

Stages 2-a (sensor and head-shape alignment plots, a QC of the conversion), 4 (filtering), 6 (epochs and evoked responses), 7 (manual MEG-MRI coregistration) and 8 (FreeSurfer) are covered on the setup page and in the pipeline’s README.md.

KIT general pipeline: environment setup and running
Stage 2: KIT .con to .fif
Stage 3-a: Bad channel detection
Stage 3-b: ICA artefact decomposition
Stage 1: trigger count sanity check (single channel)
Stage 5: build MNE events from the triggers
Kit2Fiff and coregistration GUIs: manual walkthrough

FieldTrip

First download fieldtrip from here https://www.fieldtriptoolbox.org/download/ Then, install fieldtrip folowing https://www.fieldtriptoolbox.org/download/

Fieldtrip coregistration with OPM data

You will need the headshape.pos matrix containing the head shape grid and 3 fiducials: the nasion, the lpa and the RPA (those three fiducials are at the end of the .pos file)

Remind that: Using the LaserScanner generates a basic_surface.txt file (which contains the head points) and a stylus.txt file containing 8 reference points: we will need the first three points Laser Scanner protocol reminder

Use the ft_read_head_shape to read the .pos matrix and fiducials

Information on head coordinate systems https://www.fieldtriptoolbox.org/faq/coordsys/

Beamforming, source reconstruction

  • measuring the contribution of a specific sensor in a soruce signal (weighing of the sensor w.r.t a specific source)

S(r,t) = W(r)b(t)

S rouce estimate, W beamformer weights (3 columns because 3 orientations of each dipole (add location and so on can be extended)) b channel measurement

The question is finding W (the right one among the different possible ones) adding constriants will reduce the ill-posedness of the problem - forward model constraints cfrom physiological data - inverse problem constraints

spatial filtering: assume there is only one source and estimate its activity independetly from the other sources ( the estimation is the multiplication of the weight matrix with the chanenl measurements) repeated over the number of sources

The lead matrix L is the inverse of the weight matrix W (simple assumption) Additional constraint to take the neighborhood sources into account, the W at a source r should not give high gain for another source q at a certain distance from r ( q should not be the direct neighbor) so for the direct neighbor, we wana minimise by a bit the gain and not ocmpletely set it to 0 the variance over the position is therefore minimised (ensuring smoothness, and reduced gain in neighbors) (beamformer Linearly constrained minimum variance (LCMV))

W and S are unknowns, minimising var(S) requires knowing W, but there is a relationship between W, L and C (the measured data covariance matrix, how dep/independent each measurement is with the others)

(In fieldtrip, the configuration points to LCMV or DICS(for frequency analysis) algorithm to use this beamformer)

Computing the inverse of C (used for the solution) requires that C is not rank deficient (a measure of the dependency between the rows/columns), small time window, ICA can increase deficiency

Continuous stimulus

The use of Dissimilarity matrix (is a comparison measure between two different stimulus), the matrix is low in values when the sitmulus are similar

at the stimulus level but also the brain response, computed at each time point is a way of measuring brain activity in time-continuous stimuli (video, speech,

BESA Software

The following steps are primary to process MEG data using the BESA MRI and BESA Research suite

BESA installation

Download BESA from https://www.besa.de/

The BESA license available at NYUAD-MEG lab will be soon hosted on a server, and instructions to use it will be made shortly available on this page.

You have MRI data of your participant

Open BESA MRI, start a new segmentation project, check all the segmentation options (especially BEM and FEM), pick the landmarks for segmentation and start the process. Once done, BESA will save the segmentation, BEM, FEM model outputs.

In BESA MRI, start a new coregistration project.

Open BESA Research, load your MEG data from a .fif format.

AC choosing: start with sagittal, pick the last end of the white matter fibers then go to coronal and pick the point where the white matter fibers are linked from left to right hemisphere(above the black point)

Set the nasion on the surface of the skin, because the digitization is performed on the top of the ksin

Connectivity analysis: SOme regions can be activated at the same time However, in other situations, one region fires and then triggers other regions.

It is an unsolved question to identify which begins the activity

Source montage (add to glossary): is definition of the sources

Epilepsy:can be one region activated or many other regions activating at the same time another kind of epilepsy where all the brain gets activated another kind, right hemisphere then left hemisphere

MNe kit2fiff to get a fif from .con and .mrk then open the .fif in BESA

To quickly inspect the data, and have a view of the field at a specific time, on a template head, within the sensor layout loaded from the .fif you can enable in Options -> Mapping -> Enable direct mapping by double clicking Then double click at a specific time and see the field on the template head

Radial activity within the brain cannot be picked up except for partially by the SQUIDS There is examples folder in BESA, with different datasets

Algorithms for source localization sLoreta, Loreta (named Clara in BESA, which is recursively applied), dspm Equal end current dipole (good when the phenomona is triggered by one area of the brain and not multiple areas at once) EQUILAVENT Current Dipole Beamformer (

Issues to resolve: - number of digitized head points, and alignment during coregistration - triggers in oddball experiment how to pick them up by Besa software

Open file and clear DB(to erase previous settings)

  • Identify the metrics used in computing the artifact rejection heatmap (see if we can have different metrics of these)

In artifact rejection: Low signal is to cut flat signals (that are saturated with noise)

Question about PCA components time axis

Question: Save solution works for dipoles only but not for CLARA Question: Saving activity coregistered with MRI-T1

If source analysis crashes using CLARA, reduce the voxel size

Source localization steps:

after coregistration, make sure in BESA Research, File, Coregistration, the coreg file is checked define your trials and average them, then pick a window, right click -> source analysis Pick BEM or FEM as conducting model. (the spherical model is only a template and is not used for source localization) choose a window on the left top window then pick Clara or another algorithm

Check how the USB license can be hosted on a server

Agenda: - global presentation on besa - beamformer - OPM trigger processing and coregistration

We can find the coregistration files in C:UsersPublicDocumentsBESA MRIProjectsAS-oddball-task_19010101

After source localization, you can use up and down arrows to change the sensitivity of the heatmap

SESAME uses bayesian statistics to estimate the number of sources, cortical loreta

Cortical Loreta = Loreta not on MRI but the reconstructed brain model (inflated or not inflated)

Batch-processing is the automation workflow in BESA to run similar steps multiple times

DONETODO:MATLAB import from BESA DONETODO: seeding dipole

Clara operation: goes through each time point and does the inverse (for each time point) then show the average of the time region of activation

When to use Clara or Beamformer or SESAME? Clinically: Clara (distributed model, each voxel have an activity) and Dipole fitting are approved clinically For MEG: beamformer (coz MEG has better spatial separation than EEG) SESAME (uses a dipole model, and does not assume that each voxel has an activity): DONETODO: Clara followed by SESAME, CLARA can show you regions of activity, and then SESAME can use dipoles on these regions as a prior to its operation a CLARA followed by SESAME should give a more accurate pointy result, you need to click Weight By Image (to use output of clara as input to sesame [prior])

PCA: The PCA can indicate how many sources you need, if you have 3 high activity components then you need 3 dipoles

right click a PCA component, and add to solution, this will show you where the dipole is located for that component

In Source Analysis: Residual to see what data is not covered, you can uncheck the data and keep the residual and then fit again just for the residual part

There is something called confidence level to see how the dipole explains the data (but this is nnot a validation)

In source analysis, never forget to set the baseline properly on areas where there is not much activity, prior to the stimulus

TODO: Frequency analysis difference eyes close and open, in sensor space and in source space

TODO: OPM Coregistration, how are the pink points and the sensors connected

Export NII with acivity source localized solution 1: with MIND it is a solution solution 2: .vmr file, BrainVoyager, neuroelf (in matlab) in Neuro ELFis free to import VMR, import the MRI from a .VMR file, then import analys to VMR

solution3: export after source analysis choose ACPC.nii (this setting only appears when the coregistration is set)

Solution4: longer term solution, find the transformation parameters in the project file and use them to get to ACPC coordinates, apply it on the dicom. then the exported.nii (In ACPC) activity image should match wh

OPM trigger solution:

Show code amplitude value

TODO: Send FIeldline a question about the fiducials in the .fif that has been automatically added without digitized head

Beamformer not working in oddball task because the noise level is high, (the artifacts is ok) The result can be better choosing a baseline with lower noise (-300 -200)

REgularization parameter for all methods is very important, the higher the regulalirzation parameter In beamerformer oddball, set regul parameter to 0.01 (best value), parameter accessed from Image Settings

Frequency-time analysis is not possible today with sloreta, clara, a workaround is to apply HPF and LPF filters o the region of interest then apply time-lock analysis with clara/loreta

Agenda for today:

  • finish resting state

  • do frequency analysis on time series of oddball after CLARA

==> This is not directly possible, because the orientations of the sources can be very different, in this case, the oscillation effect can double the frequency power

A directory of BESA:

  • .pdg = paradigm file (triggers, conditions, groups)

  • .fsg = averaged trials (trials and averages)

TODO(Problem): Send to BESA trigger on MISC_002 one is up and one is down, both should be up

IN resting state: - after doing an FFT, you can define your own band that your looking for in Options Band Name and Width, (you must be in SRC first) - we did beamformer in the time domain, then defined a source on the maxima obtained. then we saw the estimated time series on the maxima

BATCH creation to automate a pipeline: Shift+R or process --> Batch Processing Pause: it stops after a step in oredr for tghe user to check for things or take a screenshot, then it continues

Manual labelling of “bad” channels

Denoising

Awareness of the many sources of noise:

  • Related to the site in which the MEG system is installed

  • Related to conditions that could happen from time to time (parking garage nearby,)

Once the reasons are understood, we can identify the pattern that the noise makes.

With training data of the different possible noises, it is very possible to train a neural classifier that could identify the noise coming from the different sources and be able to denoise it from the MEG data.

Independent component analysis

Independent component analysis (ICA) is commonly used to generate what is supposed a set of independent signals from a given set of assumingly correlated signals.

The signals produced by MEG are highly correlated, therefore ICA is suitable to reduce correlation. Given a set of MEG signal X(t), ICA learns a matrix W and the output signals S(t) such that

add latex here: X(t) = W.S(t)

ICA can perform well to identify the noise signals that has a certain long lasting continuous-time pattern, but less efficient when the noise is a single event, happening at irregular periods of time.

Calling ICA withint a Python pipeline
 projs, raw.info['projs'] = raw.info['projs'], []
 ica.fit(raw)
 raw.info['projs'] = projs

Frequency Analysis

Fast-oscillating signals means high frequencies, while slow oscillations are low frequencies. In fourier space (signal represented by its Fourier transform) we can see the frequency components constituting the signal. FFT (Fast Fourier Transform) algorithm is commonly to identify the frequency components.

Research showed that signals at different frequencies have different functions at different locations of the brain. In other words, given a region of the brain, signals of frequency 8Hz are responsible of an activity that is much different than signals with frequency 20 Hz

Brain Source Estimate

When neurons become active, they do so in large groups.