Time frequency and unconstrained sources

Hello,
I have been trying to figure out which steps are the right ones when using unconstrained inversion matrices and dealing with data in the frequency domain. I have tried many ways based on tutorials and discussions here but nothing seems to make sense in my data so I thought I would ask here.

I calculated individual inversion matrices for each of my 40 participants ( MEG CTF 272 channels) based on all the epoched trials within each block (8 blocks so 8 inversion matrices per participant), I computed constrained and unconstrained but would like to use unconstrained as this was recommended when we don't have individual anatomies. In order to verify my reconstruction, I thought to check the motor lateralisation before response.

When doing so with MEG signals, I first multiplied the signals by the inversion matrice separately in the conditions where participants respond with the right hand and those where they respond with the left hand 200ms before response until respone (mean on time). This was done separately for each bloc (with the corresponding inversion matrices) then I took the norm within each bloc to go back to 15002 vertices. Then for each subject did the mean across blocs of this obtained signal. Finally I do a ttest between the condition Right/Left Hand and get a coherent lateralization signal in the sources, it is most beautiful when I do the ttest on the zscore (within each participant to account for individual variability)
image

Now I wanted to do the same in the frequency domain 8-32Hz (Donner et al. 2009 ) , I have done it before with my EEG data (El Zein et al 2015) - so I computed the fourier in these frequencies, With constrained matrices, I multiplied by the inversion matrice for all timepoints and frequencies, then took the power of abs(data) for each condition (Right/Left), then did the mean on time and frequencies, then the mean across blocs for each participant - then like described before I do a ttest of the 2 conditions and observe the result in the source space which made sense -
image

BUT I am unable to get a coherent result the unconstrained sources... I am not sure when I should take the sum of powers and what I am doing wrong. In the tutorial, it says to take the sum of the powers for each orientation, so after multiplying with the matrice and taking the abs, this is what I did - then I do the mean per block and therefore get one map per participant for each condition like previously, but there must be something wrong as the results are completely incoherent!

Could you maybe help orient the order in which I am supposed to do things with the unconstrained and frequencies?
I hope this was clear and please let me know if I need to clarify further!
Thank you very much for your help,
Marwa

The mean over time of a 200ms-segment of source signals?
This is a long time window for MEG. For example, if you had a pure sinusoid at 5Hz in these signals, it would completely disappear when averaging in time. Averaging MEG signals in time domain is not recommended over "long" periods.
Averaging the power in a specific frequency band would probably be more optimal.

(I wrote this before seeing that this was the topic of your next paragraph...)

so I computed the fourier in these frequencies, With constrained matrices, I multiplied by the inversion matrice for all timepoints and frequencies, then took the power of abs(data) for each condition (Right/Left), then did the mean on time and frequencies, then the mean across blocs for each participant - then like described before I do a ttest of the 2 conditions

Did you do all this with custom Matlab code, or did you everything with Brainstorm processes?

BUT I am unable to get a coherent result the unconstrained sources... I am not sure when I should take the sum of powers and what I am doing wrong. In the tutorial, it says to take the sum of the powers for each orientation, so after multiplying with the matrice and taking the abs, this is what I did - then I do the mean per block and therefore get one map per participant for each condition like previously, but there must be something wrong as the results are completely incoherent!

From this, I guess you tried to do it manually?
The process "Frequency > Time-frequency" and "Hilbert transform" does the TF decomposition of each signal along the 3 orientations individually, then sums the power for the 3 orientations.

For this large frequency band (8-32Hz), it is probably more indicated to use a band-pass filter + Hilbert transform, instead of this FFT+averaging in time + averaging in frequency.
Run the process "Hilbert transform" on the unconstrained sources. You may even use only one time point instead of averaging over a time window.
Does it give very different results than applying the same Hilbert transform on the constrained sources?

Thank you very much for your very fast reply!

The mean over time of a 200ms-segment of source signals?

Yes, thank you for pointing out to this, what is considered a short time window? 100ms ? 50ms?

Indeed I do everything in matlab code. All my preprocessing was done in fieldtrip then I entered the cleaned epoched files to brainstorm to compute the inversion matrices for each subject. From there I only use these imagingKernel from Brainstorm but write my own scripts, then reproject in to brainstorm to observe the effects.

The reason is that I do GLM based analysis, so I run my GLMs in the source space with my matlab scripts to get the sources for specific effects. For example at time window 200 to 300 ms after outcome the valence effect- but I was also averaging over this time course, taking the norm of the uncsontrained in the signal and then running my GLMS on each of the 150002 vertices - which may be not very accurate if I should not average over a 100ms window either?

The motor source analyses was more a way for me to confirm that my source reconstruction is valid, and when I look at the images I sent (for unconstrained MEG signals, and constrained 8-32Hz)- the lateralization does seem coherent. With the time frequency I calculated the fourier in the 8-32Hz with fieldtrip, and I use my output there that I multiply by the inversion matrices - then take the means and then take the power at each orientation for uncontrained.

For this large frequency band (8-32Hz), it is probably more indicated to use a band-pass filter + Hilbert transform, instead of this FFT+averaging in time + averaging in frequency.

By band-pass filter you mean reduce to chunks of frequencies?

Run the process "Hilbert transform" on the unconstrained sources. You may even use only one time point instead of averaging over a time window.

Do you mean after I multiply the signal by the inversion matrice (and get 45000 vertices with corresponding frequencies) I apply the Hilbert transform, before taking the power at each orientation?

Thank you very much,
Marwa

what is considered a short time window? 100ms ? 50ms?

It depends on the filters applied on your recordings, the range of frequencies you can actually observe, and what you expect in output. If you average a window of 100ms of MEG recordings, you may lose all the components above 10Hz. Computing the power in a frequency band is safer.

If you're only interested about the large transitory ERF response, it's OK. Maybe in that case you could instead low-pass filter your recordings at 10Hz, and pick only the time point in the middle of the time window you are interested in.

Indeed I do everything in matlab code.

Then unfortunately, it will be nearly impossible for us to validate what you've been doing, to reproduce erroneous behaviors, or propose improvements.

For example at time window 200 to 300 ms after outcome the valence effect- but I was also averaging over this time course, taking the norm of the uncsontrained in the signal and then running my GLMS on each of the 150002 vertices - which may be not very accurate if I should not average over a 100ms window either?

Same here. If you were low-pass filtering your signals, you would need to question how you average your time windows.

With the time frequency I calculated the fourier in the 8-32Hz with fieldtrip, and I use my output there that I multiply by the inversion matrices - then take the means and then take the power at each orientation for uncontrained

The order in which this is done in Brainstorm is tricky.
If you want to reproduce what we do without using the interface or the database, you'll have to refer directly to the code, and probably compare step-by-step the execution of your code with a real execution of bst_timefreq using the Matlab debugger:
brainstorm3/toolbox/timefreq/bst_timefreq.m at master · brainstorm-tools/brainstorm3 · GitHub

By band-pass filter you mean reduce to chunks of frequencies?

https://neuroimage.usc.edu/brainstorm/Tutorials/TimeFrequency#Hilbert_transform

Do you mean after I multiply the signal by the inversion matrice (and get 45000 vertices with corresponding frequencies) I apply the Hilbert transform, before taking the power at each orientation?

This is optimized in the code of bst_timefreq: bandpass and hilbert are linear, we can permute them with the inverse operator. We compute the complex TF coefficients for the MEG recordings, multiply them with the ImagingKernel matrix, and only then take the norm.
This is all laid out in the code of bst_timefreq.m, please refer to this function.