MVPA stc plot data extraction

If you have a question or issue with MNE-Python, please include the following info:

  • MNE-Python version: 0.23.0
  • operating system: MS Window 10

Hi everyone,

I have two questions as follows:

  1. I would like to know how I can extract the data that MNE uses to project sensor-space patterns to source space during an MVPA analysis.

Here is the code:

brain = stc.plot(hemi='split', views=('lat', 'med'), initial_time=0.1,
                 subjects_dir=subjects_dir)

and the only output that I see in my Spyder console is this:

Using control points [2.03727704 2.4722366 8.19426193]

I do not see any data from the 3D plot in my variable explorer window in Spyder. I would like to extract some data from this source model.

  1. I would also like to know what those red and blue colors mean in the source space model?

Thank you!

Best,
Aqil

:point_right: :point_right: :point_right: Please edit or remove the above text before submitting your posting. :point_left: :point_left: :point_left:

Hello @aqil.izadysadr and welcome to the forum!

The stc variable you have is an instance of SourceEstimate, and as such should contain all you need:

  • stc.time ā†’ the time points
  • stc.vertices ā†’ the vertices (left & right hemisphere)
  • stc.data ā†’ a numpy array of the data, one element per dipole and time point

Have a look at the color scale bar that should be visible on the left side of your STC plot:

red ā†’ values greater than 0
blue ā†’ values less than 0

Best wishes,

Richard

1 Like

Hi @richard,

Thank you for your response. I love your MNE tutorials on YouTube! Great Job! :+1:

And sorry for my late reply! I did not receive any notification regarding your response to my question, and I only got it today.

Regarding the second question that I asked earlier, yes the red shows values greater than 0, and blue shows values less than 0 in the source modeled MVPA. However, it looks confusing because:

  1. When you do source localization in dSPM (see Source localization with MNE/dSPM/sLORETA/eLORETA ā€” MNE 0.23.0 documentation), you do not get any of the negative values, and the source model only shows the positive ones, but why does it behave differently in source modeling MVPA and shows both the negatives and the positives?

and:

  1. Does the blue in MVPA source localization show that the activation is less than the computed mean? How about the red (activation above the computed mean)?

  2. In the following link:
    Decoding (MVPA) ā€” MNE 0.23.0 documentation
    You will see that the blues do not have negative values in the temporal generalization matrix, but they do have negative values in the source modeled MVPA, but Why?

  3. Also, what exactly does this temporal generalization matrix show in the link above?

Thank you once again!

Looking forward to hearing from you!

Best,
Aqil

The difference here is due to different parameters used when creating the inverse operator.

  • In the Source localization with MNE/dSPM/sLORETA/eLORETA example, the following parameters are used:
    inverse_operator = make_inverse_operator(
        evoked.info, fwd, noise_cov, loose=0.2, depth=0.8
    )
    
  • In the Decoding (MVPA) tutorial, on the other hand, the following parameters are used:
    inv = mne.minimum_norm.make_inverse_operator(
        evoked_time_gen.info, fwd, cov, loose=0.
    )
    

A value of loose=0 causes another parameter to be set automatically ā€“ fixed.

Letā€™s see what this fixed parameter actually does:

# %%

from pathlib import Path
import matplotlib.pyplot as plt

import mne
from mne.datasets import sample
from mne.minimum_norm import make_inverse_operator, apply_inverse

data_path = Path(sample.data_path())
raw_fname = data_path / 'MEG' / 'sample' / 'sample_audvis_filt-0-40_raw.fif'

raw = mne.io.read_raw_fif(raw_fname)  # already has an average reference
events = mne.find_events(raw, stim_channel='STI 014')

event_id = dict(aud_l=1)
tmin = -0.2
tmax = 0.5
raw.info['bads'] = ['MEG 2443', 'EEG 053']
baseline = (None, 0)
reject = dict(grad=4000e-13, mag=4e-12, eog=150e-6)

epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
                    picks=('meg', 'eog'), baseline=baseline, reject=reject)

noise_cov = mne.compute_covariance(epochs, tmax=0., method='shrunk',
                                   rank=None, verbose=True)

evoked = epochs.average()

fname_fwd = data_path / 'MEG' / 'sample' / 'sample_audvis-meg-oct-6-fwd.fif'
fwd = mne.read_forward_solution(fname_fwd)

# %%

method = 'dSPM'
snr = 3.
lambda2 = 1. / snr ** 2


def plot_stc(ax, stc, fixed):
    ax.plot(1e3 * stc.times, stc.data[::100, :].T)
    ax.set(xlabel='time (ms)', ylabel=f'{method} value')
    ax.set_title(f'fixed: {fixed}', fontweight='bold')
    

fig, axes = plt.subplots(2, 1, figsize=(6, 6), sharex=True, sharey=True)

for ax, fixed in zip(axes, (True, False)):
    inverse_operator = make_inverse_operator(evoked.info, fwd, noise_cov,
                                             fixed=fixed)
    stc = apply_inverse(evoked=evoked, inverse_operator=inverse_operator,
                        lambda2=lambda2, method=method, pick_ori=None)
    plot_stc(ax=ax, stc=stc, fixed=fixed)

axes[0].set_xlabel(None)
fig.tight_layout()

This generates the following figure:
output

For more information on this topic, I highly suggest you read through our The role of dipole orientations in distributed source localization tutorial.

Now donā€™t get confused here ā€“ these figures show very different things, and therefore similar colors donā€™t have similar meaning. Generalization across time (GAT) trains the classifier at one time point and tests it against all other time points. The resulting classification performance is shown in the GAT heatmap ā€“ blue meaning worse, red meaning better performance.All of this is done in sensor space entirely in the above example.

If you have more questions about GAT, Iā€™d suggest you open a separate topic.

Best wishes,

Richard

1 Like

Hi @richard,

Thank you for your response. Great information! I learned a lot!

In ā€œThe role of dipole orientations in distributed source localization ā€” MNE 1.7.0.dev15+gb107d92ec documentationā€, regarding the stc plot, it is mentioned:

In the plot, blue areas indicate current flowing inwards and red areas indicate current flowing outwards

So:

  1. Do inward or outward flowing currents show a different type of activation in the brain? (I am asking this to understand why we need to see them in different colors, i.e. blue and red, if they just refer to some general activation in the brain).

  2. What does that GFP, or sometimes RMS (the dotted line), which usually starts from 0, show in the MVPA source space legend? (see pictures below)

Source space with RMS:


Source space with GFP:

Thank you so much!

Best,
Aqil

Hello,

Yes and no. If by activation you mean, whether a particular region of the brain ā€œis doing somethingā€, then the answer is: it doesnā€™t matter if itā€™s blue or red; all thatā€™s important is the magnitude of the activation.

If, on the other hand, you mean, could different processes be happening in a given region if it changes color, then the answer is yes, as the overall current flow is now reversed. But keep in mind that thousands to millions of neurons might be involved here so I donā€™t think we can make any clear statements as to what exactly has changed. (But please, anyone, correct me if Iā€™m mistaken!)

It really depends on your research question and analysis approach. Keeping the signs might give you more statistical power; however, using a fixed dipole orientation (i.e. strictly orthogonal to the cortex surface) comes with drawbacks too, as described in the tutorial I shared above.

If you generate an inverse operator with ā€œlooseā€ dipole orientation and pass pick_ori='vector' to apply_inverse(), you will get a VectorSourceEstimate that still contains the 3 directional coordinates of each dipole.

If you just want to know which brain area was active and when, I would suggest to start only with the magnitudes of activation first and ignore dipole orientations, just to keep complexity low. I think it will also compute faster that way, but I might be mistaken.

RMS is the root mean square of the signal in sensor space and is essentially a measure of the ā€œagreementā€ across channels: if, for example, after a stimulus presentation the frontal channels pick up a signal thatā€™s different from what the occipital channels pick up, the RMS value will be larger than when all channels pick up similar signals. Peaks in the RMS trace are often used to find time points of interest while analyzing the data. MNEā€™s plot_joint() function, for example, finds those peaks and plots the corresponding topomaps.

GFP is global field power and calculated as the population standard deviation across all channels. Itā€™s basically the same as RMS, but with an average-reference applied to the signal first.

Because average-referencing doesnā€™t make sense in MEG, GFP is only plotted for EEG data, while RMS is plotted for MEG.

We currently donā€™t do this consistently everywhere: there might still be places where we claim to plot ā€œGFPā€ but actually show RMS; if you come across anything like this, please let us know.

Best,

Richard

1 Like

@richard

Hi Richard,

Thank you for your great response. It was very much helpful as always!

Speaking of reporting issues with MNE (while I am not quite sure if I should report the issues here), I have just noticed that when you change the time_unit parameter for the source estimate object to ā€˜msā€™ instead of the default ā€˜sā€™, the resulting source model animation does not play. The timer will keep counting, but the time bar does not move inside the Activation (AU) box.

Another thing that came to my mind is about generating videos from source estimate objects. When the source estimate object is plotted, the generate a video option does not save any videos.

Thank you!

Best,
Aqil

This should work, and if it doesnā€™t it is a bug. Please open a new bug report issue about this, including a complete code sample to load data and plot the data in a way that shows the error. Use one of the built-in datasets if possible. Please also include the output of mne sys_info (terminal command line) or import mne; mne.sys_info() (from within python console).

This also should work, and if it doesnā€™t it is a bug. Please open a separate bug report about this.

1 Like

@drammock

Hi Dan,

Thank you! Now I know how to report bugs.

Best,
Aqil