Plotting Volumetric (SEEG) activity without BEM files

[adam2392]

Describe the new feature or enhancement

It would be nice to plot volumetric data without using BEM and just the T1.mgz, or the atlas annotations from FreeSurfer.

According to #8402: "You only need a BEM if you want to do some forward modeling."

@larsoner suggested the following:

""" What you really probably want is to constrain your volumetric sources to be within the brain volume. One way to do this is by passing a BEM, but really this just uses the inner_skull surface. Instead of doing this, you can use the brain surface EDIT: by using the surf argument of setup_volume_source_space, which you can extract from the brainmask.mgz file using a few FreeSurfer commands -- this is done in the process of mne watershed_bem, so this is a convenient way to get this surface. This does not require FLASH scans -- it uses only the T1.mgz image.

Alternatively, you can tell it to use aseg.mgz or aseg+aparc.mgz with single_volume=True, which will only include volumetric sources that are in the brain volume / brainmask.mgz. """

The brainmask.mgz + T1.mgz files are obtained by running FreeSurfer recon-all for me.

Describe your proposed implementation

It seems setup_volume_source_space has an option for single_volume=True. How exactly does one use that?

If no update to API needed, then I can add some documentation to the plot_ecog/plot_seeg tutorials to make it clear how one might fit into that framework w/o the BEM files.

1. Use brain volume without BEM and demonstrate how to do this from FreeSurfer recon-all
2. update API if needed
3. update tutorial(s)

Additional comments

Stemmed from conversations in #8402

[larsoner]

It seems setup_volume_source_space has an option for single_volume=True. How exactly does one use that?

If you use single_volume=False it keeps the source spaces separate for each aseg.mgz sub-volume, but uses a lot more memory doing so. For your use case you should just use single_volume=True.

[adam2392]

Summary Workflow for SEEG (w/ and w/o BEM)

I'm going to attempt to walk you through what I've done. Please let me know if I'm just grossly misinterpreting a piece of documentation somewhere.

  # I made an evoked data structure

   # create volume source for this specific subject (not fsaverage) - I think
   fname_aseg = os.path.join(subjects_dir, subject, 'mri', 'aseg.mgz')
    vol_src = mne.setup_volume_source_space(
        subject, mri=fname_aseg,
        subjects_dir=subjects_dir,
        single_volume=True,
        verbose=True,
    )

    # create Source Time Course
    stc = mne.stc_near_sensors(
        evoked, trans, subject, subjects_dir=subjects_dir, src=vol_src,
        verbose='error')  # ignore missing electrode warnings
    clim = dict(kind='value', lims=[0, 0.4, 0.8]
                )

    # plot Nutmeg style
    fig = stc.plot(vol_src, subject=subject, subjects_dir=subjects_dir, clim=clim,
                   colormap='turbo', colorbar=True)

    # plot 3D movie and save it
    brain = stc.plot_3d(
        src=vol_src, subjects_dir=subjects_dir, colormap='turbo',
        view_layout='horizontal', views=['axial', 'coronal', 'sagittal'],
        size=(800, 300), show_traces=0.4, clim=clim,
        title=f'{bids_path.subject} {bids_path.task}',
        add_data_kwargs=dict(colorbar_kwargs=dict(label_font_size=8)))

    # save brain movie
    brain.save_movie(time_dilation=1, interpolation='linear', framerate=32,
                     time_viewer=True, filename=brain_movie_fname)

Then if one wants to use BEM files, from my understanding, you must: i) run Freesurfer recon-all , then ii) run watershed_bem using FreeSurfer, or run the corresponding mne make_watershed_bem command and then iii) feed in that subject to mne.read_source_spaces.

If one wants to instead use the brain files (w/o BEM), use setup_volume_source as done above(?)

A few issues that came up that might get answered in the meantime:

0. is the setup_volume_source correct based on what we discussed?

1. I’m a bit unclear still what is shown in the time course activation in the stc.plot_3d. I couldn’t find documentation on it in the function. Once I understand, do you think it is worth making a PR to update the doctoring?

In the plot_ecog example you stated that the time course is the It's the time course of the vertex shown by the ball, in this case the time course of the gamma activation. GFP is global field power. You can try running the example locally to play with the options and clicking around, playing with spacebar, etc.

2. Is it possible to show a time course say for a specific channel only, or specific brain region only?

3. I am currently generating a brain video but it seems that the heat map plot on the brain bleeds into the same regions. What are some tips for say constraining the source space visualization to only show the channels nearest to you?

[larsoner]

If one wants to instead use the brain files (w/o BEM), use setup_volume_source as done above(?)

Yes this should be enough to have things work

is the setup_volume_source correct based on what we discussed?

I think so, you can try swapping it in for any BEM that is used in the sEEG or ECoG example. Basically the example output should look nearly identical (maybe an extra voxel or two near the edges in the non-BEM case)

I’m a bit unclear still what is shown in the time course activation in the stc.plot_3d

It shows stc.data[idx] where idx is the index of the source with the largest absolute value across any time point. So whether it's an amplitude in current units (Am) or some statistical value or something else (here, some measure of gamma power or whatever is in evoked.data) varies depending on what is in stc.data. The critical point is that it's just this maximal source for the volume stc.

Is it possible to show a time course say for a specific channel only, or specific brain region only?

I would use epochs.copy().pick_channels(...) for the former and stc.in_label for the latter

I am currently generating a brain video but it seems that the heat map plot on the brain bleeds into the same regions. What are some tips for say constraining the source space visualization to only show the channels nearest to you?

The default is a maximum intensity projection volumetric rendering. You could change the blending from 'mip' to 'composite' in the volume_options of https://mne.tools/dev/generated/mne.VolSourceEstimate.html#mne.VolSourceEstimate.plot_3d to see if it helps. Or maybe what you really want is https://github.com/mne-tools/mne-python/issues/8382#issuecomment-712223093. It's possible tweaking the distance to be smaller or setting mode='nearest' in stc_near_sensors might help, since these also impose some smoothness.

[adam2392]

Thanks for the in depth explanation. That was very helpful.

I would like to suggest that I can make a PR into the seeg example to i) direct users to the watershed series of functions and an explanation of how to use the actual FS MRI volume/atlas if needed and ii) describe what is being shown in the time-courses plot, since it isn't immediately obvious and then iii) perhaps only showing the visualization in a certain region of the brain?

[larsoner]

i) direct users to the watershed series of functions

I guess we can. But if you don't need a BEM then I wouldn't bother computing mne watershed_bem. It's slow, a little bit flakey as a function, and produces somewhat bad surfaces for people. I would go with the mri='aseg.mgz' and volume_labels instead.

an explanation of how to use the actual FS MRI volume/atlas if needed, ii) describe what is being shown in the time-courses plot, since it isn't immediately obvious

I would point them to tutorials/source-modeling/plot_visualize_stc.py where this is covered instead of duplicating such content:

https://mne.tools/dev/auto_tutorials/source-modeling/plot_visualize_stc.html#volume-source-estimates

and then iii) perhaps only showing the visualization in a certain region of the brain?

This seems better suited to the plot_visualize_stc example since it's generally applicable beyond sEEG. So maybe put this in plot_visualize_stc as well, and in the sEEG example mention that this is one type of viz and for options they should look to plot_visualize_stc.

I'd even get rid of the stc.plot / NutMEG style one in plot_seeg in favor of stc.plot_3d, and leave the "many types of plots" to plot_visualize_stc. IMO the more we can limit the code in the sEEG and ECoG examples to things that are mostly specific to their specific considerations the better. Then we should provide links to things that sEEG and ECoG folks will want to do that are not domain specific.

So maybe a little section at the end of the example saying something like, "for additional common analyses of interest, see the following:

• For volumetric plotting options, including limiting to a specific area of the volume , see ...
• For extracting activation within a specific FreeSurfer volume, see ... etc." WDYT?

[adam2392]

Yes this makes sense. It tightens up the plot_seeg example and points them in the right direction for other types of input volume sources.

I can make a quick PR on this.