try to simulate ERD at motor cortex like in the paper

[megazero1316]

Hello , really great work and toolbox :D

however i am having some difficulty trying to simulate motor cortex ERD same as in figure 7 in the paper ( assuming that was the goal in the paper if i understand it correctly ).

so following the explanation in the paper and using the given template i tried to do the simulation on the left motor cortex

here is a sample of my code :

clear
clc
% configuring for 100 epochs of 1000 ms at 1000 Hz
t_len=1000;
config      = struct('n', 10, 'srate', 1000, 'length', t_len);

% obtaining a 64-channel lead field from ICBM-NY
leadfield   = lf_generate_fromnyhead('montage', 'S64');

% packing source locations and activations together into components
component1 = utl_get_component_fromtemplate('motorcortex_left_mu_rest_ersp', leadfield);
component2 = utl_get_component_fromtemplate('motorcortex_left_beta_rest_ersp', leadfield);
component3 = utl_get_component_fromtemplate('motorcortex_left_mu_desynch_650l_600w_ersp', leadfield);
component4 = utl_get_component_fromtemplate('motorcortex_left_beta_desynch_600l_500w_ersp', leadfield);
signal = struct('type', 'noise', 'color', 'brown','amplitude', 5, 'amplitudeDv', .5);
sourcelocs = lf_get_source_spaced(leadfield, 64, 25);
component5 = utl_create_component(sourcelocs, signal, leadfield);

% simulating data
data        = generate_scalpdata([component1,component2,component3,component4,component5], leadfield, config);
% data        = generate_scalpdata([component1,component4,component5], leadfield, config);
% data        = generate_scalpdata([component5], leadfield, config);

% converting to EEGLAB dataset format, adding ICA weights
EEG         = utl_create_eeglabdataset(data, config, leadfield);
EEG         = utl_add_icaweights_toeeglabdataset(EEG, [component1,component2,component3,component4,component5], leadfield);

and using eeglab ti plot the ERPS for the generated signal :

figure; 
pop_newtimef( EEG, 1, 28, [0  999] , [2 0.5] ...
    , 'topovec', 28, 'elocs', EEG.chanlocs, 'chaninfo', EEG.chaninfo...
    , 'caption', 'C3', 'baseline',[0]...
    , 'freqs', [3 30], 'plotitc' , 'off', 'plotphase', 'off', 'scale', 'abs', 'ntimesout', 300, 'nfreqs', 20);

i end up with this

is this correct or the way i plot it is wrong ?

[lrkrol]

Thanks for the compliment! :D

First of all I would like to express a quick warning about using the templates. I have been debating with myself whether or not to even include such a function. On the one hand, I think the general idea is good -- on the other, there's a danger that there's not sufficient generalisability to really call these "templates". Please handle with care and make sure the templates do what you want them to do.

Having said that, your procedure is generally correct, but there's a few things to keep in mind.

• As you can see, your plot ranges from about 375 to 625 ms. The desynchronisations in the template are centred around 600 and 650 ms respectively, so they won't show up very clearly in this plot. pop_newtimef requires quite a bit of time points for its calculations, which is why so much of the data range is missing. In the paper, I simulated a 3-second epoch, thus having enough "padding" data for this function.

• The output of pop_newtimef is generally quite sensitive to its various parameters. IIRC, the config used in the paper was:

figure;
pop_newtimef( EEG, 1, 28, [-1000  1999], [2 0.1], ...
     'topovec', 28, 'elocs', EEG.chanlocs, 'chaninfo', EEG.chaninfo, ...
     'caption', 'C3', 'baseline',[0],'basenorm', 'off', 'nfreqs', 28, 'freqs', [3 30], ...
     'plotitc' , 'off', 'plotphase', 'off', 'padratio', 32, 'freqscale', 'log');

• You are adding both resting-state mu/beta and desynchronisations. The "rest" components project a steady signal with power in the mu/beta band. The desynchronisation components also project that same signal, except it is reduced to half-amplitude during the desynch window. Adding both rest and desynch components (on the same side) thus increases the power overall and makes the desynch less significant. You probably wanted to put the rest components on the right side, not the left, or vice versa for the desynchs.

The following adjusted code increases the length of the epoch (also adding a 1-second pre-stimulus interval) and takes only components 3, 4, and 5. Then using the above code should plot a figure making the desynch more clearly visible.

% configuring for 100 epochs of 1000 ms at 1000 Hz
t_len=3000;
prestim = 1000;
config      = struct('n', 10, 'srate', 1000, 'length', t_len, 'prestim', prestim);

% obtaining a 64-channel lead field from ICBM-NY
leadfield   = lf_generate_fromnyhead('montage', 'S64');

% packing source locations and activations together into components
component3 = utl_get_component_fromtemplate('motorcortex_left_mu_desynch_650l_600w_ersp', leadfield);
component4 = utl_get_component_fromtemplate('motorcortex_left_beta_desynch_600l_500w_ersp', leadfield);
signal = struct('type', 'noise', 'color', 'brown','amplitude', 5, 'amplitudeDv', .5);
sourcelocs = lf_get_source_spaced(leadfield, 64, 25);
component5 = utl_create_component(sourcelocs, signal, leadfield);

% applying prestimulus
component3 = utl_shift_latency(component3, prestim);
component4 = utl_shift_latency(component4, prestim);

% simulating data
data        = generate_scalpdata([component3,component4,component5], leadfield, config);

% converting to EEGLAB dataset format, adding ICA weights
EEG         = utl_create_eeglabdataset(data, config, leadfield);
EEG         = utl_add_icaweights_toeeglabdataset(EEG, [component3,component4,component5], leadfield);

(When writing this, I noticed a bug in utl_shift_latency; please re-download or fix locally!)

Some more things to keep in mind:

• The plot in the paper was from an ICA cluster, not from raw scalp data. The ICA cluster really only showed the motor cortex component activity, no noise components. Scalp data, which you are plotting, of course has both.

• 10 epochs is generally a very small number to get reliable effects. Increasing this number will further average out the noise, making the desynch more obvious. (Also see the section in the tutorial about adjusting signal to noise ratios.)

When increasing the number of epochs to 50 using the above code, I get the following plot:

[megazero1316]

Thanks for your help , now i see where is my mistakes 😄 .

the reason why i plot raw scalp data is that i train some machine learning models on bci competition iv dataset 2a which it`s for motor imagery EEG signal. and i plan to simulate some EEG signal using your toolbox to validate that these models focus on picking up ERD/ERS features, so what i do is i visualize the ERD/ERS for a subject and i want to try to simulate similar pattern using your toolbox.

so i have a question, is there anything i should consider when i do the simulation to achieve a toy data similar to one of the subject (e.g. source location and dipole orientation variation, the deviation in frequency and latency and so on) ?

[lrkrol]

I've never really worked with motor imagery datasets myself, but I think participants usually sustain the imagery for a longer amount of time (e.g. multiple seconds). The one that I simulated in the paper, and thus the one that's in the template file, is a rather short ERD that mimicks an actual manual response (e.g. a single button press).

So, you may want to simulate longer epochs with a more sustained ERD/ERS effect. The template uses the invburst ERSP modulation to get an ERD at a specific time point. In a simple motor imagery simulation for classification, you may not even need that -- you can just simulate one set of epochs that have a relatively strong signal in the mu band coming from some relevant motor cortex dipole, and another set of epochs with lower power in that band. You could use frequencyShift to have the signal amplitudes between classes be different on average but overlap so as not to make it too easy.

That would probably be the simplest case, good for classification, but less so for visualisation, because it wouldn't have a baseline period. If you want to have a baseline and "ramp-up" period, the invburst is probably fine, but as I said, your real data probably has a longer-lasting effect, so you'd want to increase the width of the burst.

If you want to make toy data similar to one of your real datasets, I would suggest to look at that real data and get the information you need. Perhaps you could obtain some single-trial statistics to see the mean and standard deviation (e.g. of the signal power) in the classes and use those to inform the simulation parameters. As for the location of the source, you could have a look at the projection patterns in the mu band (or CSP patterns, if you're using CSP) and try to reproduce them visually, or apply a dipole-fitting method.

Also note that the variation in dipole orientation may not be realistic. I included it because it's technically possible, but anatomically, the cortex obviously doesn't do a lot of rotating on the spot. I would initially stick to amplitude and perhaps latency variations in this case, and make sure there is a realistic amount of noise (also see the SNR section of the readme).