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
Follow instructions here MNE Install Ideally, choose the standalone installer it usually has the complete suite. We recommend to install MNE as a standalone installer.
The data of the MNE environment will be under C:\ProgramData\mne-python\1.x.y_z (for Windows)
You can setup Pycharm (or your favorite IDE) to use this environment as a Python environment for your pipeline project.
The configuration file for setting SUBJECTS_DIR, the directory to the data of your subject, can be set in:
C:\Users\user_name\.mne\mne-python.json (for Windows)
Installing freesurfer on windows (For source localization with MNE)
If you want to do source localization in MNE, you will need to install freesurfer.
If you are under windows, configure WSL2 on windows. You can then access the files of the ubuntu distribution by typing \\wsl$\Ubuntu in the file explorer in windows.
Install freesurfer following their documentation page. https://surfer.nmr.mgh.harvard.edu/fswiki/FS7_wsl_ubuntu
Order of notebooks for general MNE-KIT pipelines
You have acquired MEG-KIT data from NYUAD and wish to start your pre-processing of your data:
If your experiment design is based on using single-channel mode` for coding your triggers, use the following sequence of pipelines to:
convert the raw KIT .con files (with and without CALM noise reduction processing) into .fif
[Optional but highly recommended] perform trigger count sanity check (trigger count and trigger sequence check), this pipeline ensures no triggers are lost and that they are detected correctly
compute event sample onset and save .eve files
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.
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.