Example analysis pipelines

The optimal methods for cleaning TMS-EEG data are an active areas of research and, as such, are continually evolving. Here, we give one example of a possible pipeline for TMS-EEG analysis. This can be considered a perpetual beta version and comes with no guarantees. The basic steps for analysis are summarised below. We then describe how to complete each step from the EEGLAB user interface and how to output this pipeline to a script using the eegh command. For further detail on the theoretical rationale behind these steps, please see the following paper:

• Paper coming soon.

Example dataset and script

An example dataset at several different stages of processing is available to test this pipeline and can be downloaded from here: https://figshare.com/articles/TESA_example_data_and_scripts/3188800

An example script demonstrating the following pipeline is available within the TESA toolbox: example_script_from_manual.m

Analysis pipeline steps.

1. Import data to EEGLAB

2. Load channel locations

3. Remove unused electrodes

4. Find TMS pulses (if no triggers recorded)

5. Remove bad electrodes

6. Epoch data (-1000 to 1000 ms)

7. Demean data (-1000 to 1000 ms) ‡

8. Remove TMS pulse artifact and peaks of TMS-evoked muscle activity (-2 to 10 ms)

9. Interpolate missing data around TMS pulse ^

10. Downsample data (5000 Hz to 1000 Hz)

11. Remove bad trials

12. Replace interpolated data around TMS pulse with constant amplitude data (-2 to 10 ms) † ~

13. Remove large amplitude artifacts including TMS-evoked muscle, electrical, and movement artifacts (using FastICA and auto component selection) †

14. Extend data removal to 15 ms (-2 to 15 ms) †

15. Interpolate missing data around TMS pulse † ^

16. Bandpass (1-100 Hz) and bandstop (48-52 Hz) filter data

17. Replace interpolated data around TMS pulse with constant amplitude data ~

18. Remove all other artifacts (using FastICA an auto component selection)

19. Interpolate constant amplitude data around TMS pulse

20. Interpolate missing channels

21. Extract a region of interest

22. Find peaks of interest

23. Output peak analysis in table

24. Plot the results

‡ Baseline correction can reduce ICA reliability. Demeaning the data is an alternative at this step to remove DC-offsets.

† These steps are only necessary if TMS-evoked muscle, electrical, or movement artifacts are present in the data.

^ Interpolating missing data following the removal of the TMS pulse is necessary prior to filtering to avoid sharp steps or transitions in the data, which can lead to ringing artifacts.

~ Replacing interpolated data around TMS pulse with constant amplitude data is necessary prior to ICA to improve performance.

Step 1: Import data to EEGLAB

Raw data from the acquisition software and EEGLAB data are imported in two different ways:

From raw data file (.vhdr):

File -> Import data -> Using EEGLAB functions and plugins -> Brain Vis. Rec. .vhdr file -> select file -> OK.

From EEGLAB data file (.set - example data set provided):

File -> Load existing dataset -> select file.

Step 2: Load channel locations

Load channel locations for plotting topographic maps. This is required by several EEGLAB and TESA functions.

Edit -> Channel locations -> Press OK for standard positions (if you have recorded subject specific locations, these can also be loaded here) -> Press OK.

Step 3: Remove unused electrodes

If electrodes were recorded from but are not used, remove them here.

Edit -> Select data -> next to 'channel range' click the tick box -> press the grey button across from channel range -> hold ctrl and select electrodes 31 and 32 (unused electrodes) -> press OK -> press OK.

Step 4: Find TMS pulses (if no triggers recorded)

Note there are triggers in this file already so can skip this step.

Tools -> TMS-EEG signal analyser (TESA) -> Find TMS pulse -> change refractory time to 4 (may not be required for other data) -> Press OK.

Step 5: Remove bad electrodes

Here we make use of one of EEGLAB's functions that automatically detects bad electrodes based on kurtosis.

Tools -> Automatic channel rejection -> Press OK -> When plot appears check electrodes and press reject -> Press OK.

Step 6: Epoch data (-1000 to 1000 ms)

Wide epochs are recommended to provide adequate time from the period of interest (0-500 ms) to prevent boundary artifacts when filtering.

Tools -> Extract epochs -> Press the grey box next to 'time-locking event type' (use either R128 or TMS) -> Change 'Epoch times' to -1 1 -> Press OK

Step 7: Demean data (-1000 ms to 1000 ms)

Use the remove baseline function to subtract the average of the entire epoch from each data point. Note the baseline correction window will automatically generate after epoch data. If cancel is pressed, the window can be accessed from the user interface via the following:

Tools -> Remove baseline -> Change 'baseline latency range' to -1000 1000 -> Press OK.

Step 8: Remove TMS pulse artifact and peaks of TMS-evoked muscle activity (-2 to 10 ms)

Remove the TMS pulse artifact and replace with 0s.

Tools -> TMS-EEG signal analyser (TESA) -> Remove artifact data -> Ensure 'Time to remove TMS artifact' is -2,10 -> Press OK.

Step 9: Interpolate missing data around TMS pulse

Use cubic interpolation to replace the removed data. This helps minimise ringing or step artifacts caused by the low-pass, anti-aliasing filter applied prior to downsampling. The aim here is to create a smooth transition between the real and interpolated data. The values chosen for the 'Time before removed data' and 'Time after removed data' included for fitting the cubic function will depend on the sampling rate of your data and how 'sharp' the transition is between 0 and your real data. Note that at least two samples must be included in these windows for cubic fitting to work (e.g. 1 ms @ 1 kHz only has one sample, whereas 1 ms at 5 kHz has 5 samples, as is the case in this example data set).

Tools -> TMS-EEG signal analyser (TESA) -> Interpolate removed data -> Selec 'cubic' from the dropdown menu -> Change the 'Time before removed data' and 'Time after removed data' to 1 -> Press OK.

Step 10: Downsample data (5000 Hz to 1000 Hz)

Reduce the sampling rate to reduce the file size. A high sampling rate will significantly slow computational time and is not required as we will low pass filter to 100 Hz.

Tools -> Change sampling rate -> New sampling rate = 1000 -> Press OK.

Step 11: Remove bad trials

Use one of EEGLAB's functions to automatically detect bad trials based on the probability of the data and then manually check the results. Note that the SD thresholds can be altered here if 5 SDs is marking too many/few trials.

Tools -> Reject data epochs -> Reject data (all methods) -> Under 'Find improbable data' click calculate -> Click 'Scroll data' in the top right corner -> Manually scroll through the data and check the epochs marked for rejection. When finished, click update marks (bottom right corner) -> Press OK. -> On the 'Reject trials using data statistics' window, click Reject marked trials (bottom right). -> Click yes -> Press OK.

Step 12: Replace interpolated data around TMS pulse with constant amplitude data (-2 to 10 ms)

Again remove the interpolated data. Replacing interpolated data around TMS pulse with constant amplitude data is necessary prior to ICA to improve performance.

Tools -> TMS-EEG signal analyser (TESA) -> Remove artifact data -> Ensure 'Time to remove TMS artifact' is -2,10 -> Press OK.

Step 13: Remove large amplitude artifacts including TMS-evoked muscle, electrical, and movement artifacts (using FastICA and automatic component selection)

Use one of the approaches for removing TMS-evoked muscle artifacts. We have used FastICA with automatic component selection, however EDM, PCA compression and suppression, or detrending could be used here instead.

Tools -> TMS-EEG signal analyser (TESA) -> FastICA -> Press OK.

Wait for FastICA to finish.

Tools -> TMS-EEG signal analyser (TESA) -> Auto component selection -> Uncheck Eyeblinks, Lateral eye movement, Persistent muscle activity, Electrode noise -> Enter 15 next to 'Components to plot' -> Press OK. -> Check each component classification and press enter when satisfied. -> Check pre and post component rejection plots and press enter when satisified -> Indicate whether to continue (Yes) or repeat component selection (No- Redo).

In the example data, there is one component which comes up consistently that is not a TMS-evoked muscle artifact, but is worth removing in the first run of FastICA. This component represents a large step artifact, which could prove problematic for filtering in the next stage of analysis. Manually mark this artifact as 'other' and it will be removed.

Step 14: Extend data removal to 15 ms (-2 to 15 ms)

This is to remove any residual muscle activity not removed above.

Tools -> TMS-EEG signal analyser (TESA) -> Remove artifact data -> Ensure 'Time to remove TMS artifact' is -2,15 -> Press OK.

Step 15: Interpolate missing data around TMS pulse

Use cubic interpolation to replace the removed data. This helps minimise ringing or step artifacts caused by the bandpass/bandstop filtering. Note that because we have downsampled the data to 1 kHz, we now use a 5 ms window, which included 5 samples, for fitting the cubic function.

Tools -> TMS-EEG signal analyser (TESA) -> Interpolate removed data -> Selec 'cubic' from the dropdown menu -> Change the 'Time before removed data' and 'Time after removed data' to 5 -> Press OK.

Step 16: Bandpass (1-100 Hz) and bandstop (48-52 Hz) filter data

Filtering helps remove non-neural data and may improve ICA decomposition.

Tools -> TMS-EEG signal analyser (TESA) -> Butterworth filter -> Check frequencies are 1 and 100, filter order is 4 and filter type is band-pass -> Press OK.

Tools -> TMS-EEG signal analyser (TESA) -> Butterworth filter -> Change frequencies to 48 and 52, and filter type is band-stop -> Press OK.

Step 17: Replace interpolated data around TMS pulse with constant amplitude data (-2 to 15 ms)

Again remove the interpolated data. Replacing interpolated data around TMS pulse with constant amplitude data is necessary prior to ICA to improve performance.

Tools -> TMS-EEG signal analyser (TESA) -> Remove artifact data -> Ensure 'Time to remove TMS artifact' is -2,15 -> Press OK.

Step 18: Remove all other artifacts (using FastICA and automatic component selection)

Use FastICA and automatic component selection to remove all other types of artifacts such as blinks, eye movement, persistent muscle activity and electrode noise.

Tools -> TMS-EEG signal analyser (TESA) -> FastICA -> Press OK.

Wait for FastICA to finish.

Tools -> TMS-EEG signal analyser (TESA) -> Auto component selection -> Press OK (to use default settings). -> Check each component classification and press enter when satisfied. -> Check pre and post component rejection plots and press enter when satisified -> Indicate whether to continue or repeat component selection.

Step 19: Interpolate missing data around TMS pulse

Interpolate the removed data. At this stage this is largely aesthetic.

Tools -> TMS-EEG signal analyser (TESA) -> Interpolate removed data -> Press OK.

Step 20: Interpolate missing channels

Interpolate any missing channels that were removed at the start of the pipeline using spherical interpolation.

Tools -> Interpolate electrodes -> Select 'Use all channels from other dataset' -> In 'Dataset index' enter a stored dataset which contains all of the original electrodes prior to removal (e.g. 2). -> Press OK.

Step 21: Re-reference to average

We recommend an average reference, but this may not be possible for sparse recording montages. Note that in this pipeline we only average reference towards the end of the pipeline. This is because we have removed certain channels earlier in the pipeline, which results in an uneven distribution of channels, invalidating one of the assumptions of average referencing. As such, we only average reference once these removed channels have been interpolated. However, it is unclear whether this issue of uneven electrode distributions and average referencing really matters from a practical stand point.

Tools -> Re-reference -> Ensure 'Compute average re-reference' is selected -> Press OK.

Step 22: Extract a region of interest

Extract a region of interest which will average over trials. This could be a single electrode or the average of a group of electrodes.

Tools -> TMS-EEG signal analyser (TESA) -> Extract TEPs -> Select electrodes -> Hold Ctrl and select P1, P3, CP1 and CP3 -> Press OK -> Give the ROI the name Parietal next to 'Identifier for ROI analysis' -> press OK.

Step 23: Find peaks of interest

Look for three positive peaks at 40, 80 and 200 ms and three negative peaks at 20, 60 and 100 ms. The time windows for each of the peaks are: 20 (10-30 ms), 40 (30-50 ms), 60 (50-70 ms), 80 (70-90 ms), 100 (90-110 ms), 200 (180-220 ms).

Positive peaks

Tools -> TMS-EEG signal analyser (TESA) -> Find and analyse TEP peaks -> At 'Input type for ROI analysis' select ROI.Parietal -> At 'Direction of peaks' select positive -> At 'Peaks to find' type: 40, 80, 200 -> At 'Time window for peak search' type: 30,50; 70,90; 180,220 -> Press OK.

Negative peaks

Tools -> TMS-EEG signal analyser (TESA) -> Find and analyse TEP peaks -> At 'Input type for ROI analysis' select ROI.Parietal -> At 'Direction of peaks' select negative -> At 'Peaks to find' type: 20, 60, 100 -> At 'Time window for peak search' type: 10,30; 50,70; 90,110 -> Press OK.

Step 24: Output peak analysis in table

Return the peak analysis results in a table. This can then be copied and pasted in to other programs if required.

Tools -> TMS-EEG signal analyser (TESA) -> Output TEP peak analysis -> At 'TEP for output' select ROI.Parietal -> At 'Calculation type' select amplitude -> At 'Peaks for calculation' select indiviual -> Press OK.

Step 25: Plot the results

Plot the results of the ROI analysis.

Tools -> TMS-EEG signal analyser (TESA) -> Plot data -> At 'Data to plot' select ROI parietal -> At 'Plot peak analysis results' click on tick box -> Press OK.

Convert the manual analysis in to a script

Once the manual analysis is complete through the EEGLAB user interface, the command line functions can be generated using the EEGLAB history. These can then be pasted in to a script which can be run to semi-automate the analysis.

Use eegh to output function syntax

Type eegh in to the command line and press enter. The following will be returned in the command window. This is a list of all of the functions run by EEGLAB during the preceeding analysis. These functions can be copied and pasted in to a script, which allows the user to repeat the analysis without using the EEGLAB interface. Under 'Home' press 'New Script' to open the editor and paste these functions in to the script. This script can now be saved for future use.

[ALLEEG EEG CURRENTSET ALLCOM] = eeglab;

EEG = pop_loadset('filename','example_data.set','filepath','H:\TESA_example_data\test_TESA\');

[ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 );

EEG=pop_chanedit(EEG, 'lookup','C:\Program Files\MATLAB\eeglab13_5_4b\plugins\dipfit2.3\standard_BESA\standard-10-5-cap385.elp');

[ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET);

EEG = eeg_checkset( EEG );

EEG = pop_select( EEG,'nochannel',{'31' '32'});

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 1,'gui','off');

EEG = eeg_checkset( EEG );

EEG = pop_rejchan(EEG, 'elec',[1:62] ,'threshold',5,'norm','on','measure','kurt');

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 2,'gui','off');

EEG = eeg_checkset( EEG );

EEG = pop_epoch( EEG, { 'R128' }, [-1 1], 'epochinfo', 'yes');

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 3,'gui','off');

EEG = eeg_checkset( EEG );

EEG = pop_rmbase( EEG, [-1000 1000]);

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 4,'gui','off');

EEG = pop_tesa_removedata( EEG, [-2 10] );

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 5,'gui','off');

EEG = pop_tesa_interpdata( EEG, 'cubic', [1 1] );

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 6,'gui','off');

EEG = eeg_checkset( EEG );

EEG = pop_resample( EEG, 1000);

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 7,'gui','off');

EEG = eeg_checkset( EEG );

EEG = pop_jointprob(EEG,1,[1:59] ,5,5,0,0);

EEG = pop_tesa_removedata( EEG, [-2 10] );

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 8,'gui','off');

EEG = pop_tesa_fastica( EEG, 'approach', 'symm', 'g', 'tanh', 'stabilization', 'off' );

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 9,'gui','off');

EEG = pop_tesa_compselect( EEG,'compCheck',on,'comps',15,'figSize','small','plotTimeX',[-200 500],'plotFreqX',[1 100],'tmsMuscle','on','tmsMuscleThresh',8,'tmsMuscleWin',[11 30],'tmsMuscleFeedback','off','blink','off','blinkThresh',2.5,'blinkElecs',{'Fp1','Fp2'},'blinkFeedback','off','move','off','moveThresh',2,'moveElecs',{'F7','F8'},'moveFeedback','off','muscle','off','muscleThresh',0.6,'muscleFreqWin',[30 100],'muscleFeedback','off','elecNoise','off','elecNoiseThresh',4,'elecNoiseFeedback','off' );

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 10,'gui','off');

EEG = pop_tesa_removedata( EEG, [-2 15] );

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 11,'gui','off');

EEG = pop_tesa_interpdata( EEG, 'cubic', [5 5] );

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 12,'gui','off');

EEG = pop_tesa_filtbutter( EEG, 1, 100, 4, 'bandpass' );

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 13,'gui','off');

EEG = pop_tesa_filtbutter( EEG, 48, 52, 4, 'bandstop' );

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 14,'gui','off');

EEG = pop_tesa_removedata( EEG, [-2 15] );

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 15,'gui','off');

EEG = pop_tesa_fastica( EEG, 'approach', 'symm', 'g', 'tanh', 'stabilization', 'off' );

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 16,'gui','off'); EEG = pop_tesa_compselect( EEG,'compCheck','on','comps',[],'figSize','small','plotTimeX',[-200 500],'plotFreqX',[1 100],'tmsMuscle','on','tmsMuscleThresh',8,'tmsMuscleWin',[11 30],'tmsMuscleFeedback','off','blink','on','blinkThresh',2.5,'blinkElecs',{'Fp1','Fp2'},'blinkFeedback','off','move','on','moveThresh',2,'moveElecs',{'F7','F8'},'moveFeedback','off','muscle','on','muscleThresh',0.6,'muscleFreqWin',[30 100],'muscleFeedback','off','elecNoise','on','elecNoiseThresh',4,'elecNoiseFeedback','off' );

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 17,'gui','off');

EEG = pop_tesa_interpdata( EEG, 'cubic', [5 5] );

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 18,'gui','off');

EEG = eeg_checkset( EEG );

EEG = pop_interp(EEG, ALLEEG(2).chanlocs, 'spherical');

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 19,'gui','off');

EEG = eeg_checkset( EEG );

EEG = pop_reref( EEG, []);

[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 20,'gui','off');

EEG = pop_tesa_tepextract( EEG, 'ROI', 'elecs', {'P3','CP1','P1','CP3'}, 'tepName', 'Parietal' );

EEG = pop_tesa_peakanalysis( EEG, 'ROI', 'positive', [40 80 200], [30 50;70 90;180 220], 'method', 'largest', 'samples', 5 );

EEG = pop_tesa_peakanalysis( EEG, 'ROI', 'negative', [20 60 100], [10 30;50 70;90 110], 'method', 'largest', 'samples', 5 );

output = pop_tesa_peakoutput( EEG, 'tepName', 'Parietal', 'calcType', 'amplitude', 'winType', 'individual', 'averageWin', 5, 'fixedPeak', [], 'tablePlot', 'on' );

pop_tesa_plot( EEG, 'tepType', 'ROI', 'tepName', 'Parietal', 'xlim', [-100 500], 'ylim', [], 'CI','off','plotPeak','on' );

Alter the script to improve use.

An altered version of the above pipeline is available in the TESA toolbox: example_script_from_manual.m

To run, load the script in to the Matlab editor and press the green run button under the editor tab. The script will then sequentially run through each function.

The following modifications were made to make the script run more smoothly.

• Remove all of the following functions: eeg_store, pop_newset and eeg_checkset. These functions are not necessary when running scripts as the data does not need to be stored for later access. It is better to save them to file (see below).

• Add comments to annotate each step. Placing a % in front of a line in a script means this line is considered a comment and will not be evaluated. This is useful for annotation.

• Add clear; close all; clc; to the top of the script. This clears all data from the workspace, closes any figures that are open and clears the command window. It's always good to start with a clean slate.

• Use substitution and concatenation to set file name and file path for loading. Substitution involves saving the file path name as a string (e.g. text) variable in the workspace (coded as: filepath = 'H:\\TESA_example_data\\test_TESA\\';). Instead of entering the full path as text in the load function, the variable filepath can now be entered instead. The same can be applied for the file name (coded as: fileprefix = 'example_data';). Concatenation involves combining several strings together to make one string. This is useful for building file names (see save points below) and is coded as: filename = [fileprefix, '.set']; The filename variable returns the following string: 'example_data.set'. Now if either the path where the data is stored or the name of the data changes (e.g. if you want to apply the same script to a different data file), changing the filepath and fileprefix at the top of the script will automatically apply these new names to the rest of the script.

• Save the original channel locations in the EEG structure and use this saved data for the interpolation step. After the channel locations have been loaded (stored in the EEG structure in the workspace under EEG.chanlocs) and unused channels removed, include a line which saves the original channel locations for later use in the interpolation function (coded as: EEG.allchan = EEG.chanlocs;). In the interpolation function (pop_interp), use this variable as input for the original channel locations (replace ALLEEG(2).chanlocs with EEG.allchan).

• Generalise the remove bad channel step for different sized electrode arrays. Change the pop_rejchan to the following: EEG = pop_rejchan(EEG, 'elec',[1:size(EEG.data,1)] ,'threshold',5,'norm','on','measure','kurt');. The command size(EEG.data,1) will return the size of the first dimension of the EEG data matrix, which corresponds to the electrode number. Therefore, if this script is used on a different data set with a different number of electrodes (e.g. not 62), it will still work.

• Alter the bad trial detection to include a manual check when running with a script. Change pop_jointprob to the following: EEG = pop_jointprob(EEG,1,[1:size(EEG.data,1)] ,5,5,0,0);. This will adapt the function to look for the number of channels present in the data. Change pop_rejepoch to EEG = pop_rejepoch( EEG, EEG.BadTr ,0);. Above this line create a new line and enter the following: EEG.BadTr = unique([find(EEG.reject.rejjp==1) find(EEG.reject.rejmanual==1)]);. These two lines will now save trials that have been automatically detected and manually marked in the EEG structure as EEG.BadTr and will then reject these epochs. Next, replace EEG = eeg_rejsuperpose( EEG, 1, 1, 1, 1, 1, 1, 1, 1); with pop_rejmenu(EEG,1);. This will plot the reject trials screen. Underneath this line, create a new line as follows: pause_script = input('Highlight bad trials, update marks and then press enter');. This line will pause the script while you manually check through the trials (by pressing scroll data in the top right corner). When you have finished, press update marks and close the main trial rejection window (instead of pressing Reject marked trials). Now press enter in the command line and the script will continue. This is a bit clunky, but is necessary if you would like to manually check the trials. Note: The input function has changed in Matlab R2016a and will no longer allow the user to interactively access the figure while the script is paused. This works for Matlab R2015b and earlier.

• Include some save points. We've included two save points after cleaning and after analysis. It is good practice to regularly save at different cleaning stages, so more save points could be included as required. First use concatenation to create a new filename so that the existing data file is not saved over. filename = [fileprefix,'_cleaned.set']; Here the filename will be example_data_cleaned.set. Then include a line for saving using the EEGLAB save function. EEG = pop_saveset( EEG, 'filename',filename,'filepath',filepath);

Note that the file path to the default electrode positions in pop_chanedit may need to be altered depending on where the eeglab folder is stored.