mne-python Extracting the time series of activations in a label

Describe the new feature or enhancement

Current API: https://mne.tools/stable/generated/mne.SourceEstimate.html#mne.SourceEstimate.extract_label_time_course implementation offers 4 modes to extract source time courses in a label.

Modes: max, mean, mean_flip, pca_flip

Is it possible to extend the mode options to com and min modes?

com: time course of the vertex that represents spatial center of mass in a label.
min: minimum value across vertices at each time point within each label.

Reasons:

For dSPM source time-courses with singed values, sometimes min mode (Minimum value across vertices) is more representative than max mode if the underlying label (dipoles) represent a surface negative source; for example: auditory cortex (Heschl's gyrus: N1 component). With a prior information (if the source/label is well known: surface positive/negative), min/max mode can be more useful than any other strategies to average the times series in a label.
com: dipole that represent spatial center of mass in a label. com can be an extra option to extract source time-courses in a label i.e., shows similar values as min/max modes.

Example to follow: https://mne.tools/stable/auto_examples/inverse/label_source_activations.html#sphx-glr-auto-examples-inverse-label-source-activations-py

Current, max mode in MNE is implemented as:

_label_funcs = {
    "mean": lambda flip, data: np.mean(data, axis=0),
    "mean_flip": lambda flip, data: np.mean(flip * data, axis=0),
     "max": lambda flip, data: np.max(np.abs(data), axis=0),
    "pca_flip": _pca_flip,
}

isn't it make more sense not to take the max(absolute(data)) and rather use only max (data) especially, if I consider from the signed prospective? I wrote some code to test this.

Describe your proposed implementation

test code:

import matplotlib.pyplot as plt
import matplotlib.patheffects as path_effects

import mne
from mne.datasets import sample
from mne.minimum_norm import read_inverse_operator, apply_inverse
import numpy as np

data_path = sample.data_path()
label = "Aud-lh"
meg_path = data_path / "MEG" / "sample"
subjects_dir = data_path / "subjects"
label_fname = meg_path / "labels" / f"{label}.label"
fname_inv = meg_path / "sample_audvis-meg-oct-6-meg-inv.fif"
fname_evoked = meg_path / "sample_audvis-ave.fif"

snr = 3.0
lambda2 = 1.0 / snr**2
method = "dSPM"  # use dSPM method (could also be MNE or sLORETA)

# Load data
evoked = mne.read_evokeds(fname_evoked, condition=0, baseline=(None, 0))
inverse_operator = read_inverse_operator(fname_inv)
src = inverse_operator["src"]

# %%
# Compute inverse solution
# ------------------------
pick_ori = "normal"  # Get signed values to see the effect of sign flip
stc = apply_inverse(evoked, inverse_operator, lambda2, method, pick_ori=pick_ori)

label = mne.read_label(label_fname)

stc_label = stc.in_label(label)
modes = ("com", "mne_max", "min", "max",  "mean", "mean_flip", "pca_flip")
tcs = dict()
for mode in modes:
    if mode == "com":
        # compute spatial center of mass (COM)
        label.values.fill(1.0)

        # Restrict the eligible vertices to be those on the surface under
        # consideration and within the label.
        hemi_to_ind = {"lh": 0, "rh": 1}
        surf_vertices = src[hemi_to_ind[label.hemi]]["vertno"]
        restrict_verts = np.intersect1d(surf_vertices, label.vertices)
        com = label.center_of_mass(subjects_dir=subjects_dir, subject='sample', restrict_vertices=restrict_verts)
        com_dipole = np.where(stc_label.vertices[0] == com)[0][0]
        tcs[mode] = [stc_label.data[com_dipole, :]]

    elif mode == "mne_max":  # current implementation in mne
        tcs[mode] = [np.max(abs(stc_label.data), axis=0)]

    elif mode == "min":
        tcs[mode] = [np.min(stc_label.data, axis=0)]

    elif mode == "max":
        tcs[mode] = [np.max(stc_label.data, axis=0)]

    else:
        tcs[mode] = stc.extract_label_time_course(label, src, mode=mode)
print("Number of vertices : %d" % len(stc_label.data))

# %%
# View source activations
# -----------------------

fig, ax = plt.subplots(1)
t = 1e3 * stc_label.times
ax.plot(t, stc_label.data.T, "k", linewidth=0.5, alpha=0.5)
pe = [
    path_effects.Stroke(linewidth=5, foreground="w", alpha=0.5),
    path_effects.Normal(),
]
for mode, tc in tcs.items():
    ax.plot(t, tc[0], linewidth=3, label=str(mode), path_effects=pe)
xlim = t[[0, -1]]
ylim = [-30, 30]
ax.legend(loc="upper right")
ax.set(
    xlabel="Time (ms)",
    ylabel="Source amplitude",
    title="Activations in Label %r" % (label.name),
    xlim=xlim,
    ylim=ylim,
)
mne.viz.tight_layout()
plt.show()

Test result: label_activations

Describe possible alternatives

None

Additional context

No response

Aug 24 '23 16:08 dasdiptyajit

how about allowing the mode parameter to be a callable so you can do any logic you need?

Message ID: @.***>

Aug 25 '23 07:08 agramfort

Some thoughts

docs are wrong about max, they don't mention that it's max absolute value
changing the semantics of any of the existing options is probably a bad idea, too many people have probably memorized their preferred option and/or copy-paste old code; we don't want the behavior to surprise them / have them misinterpret the finding because they didn't know the behavior changed.
although min and (non-absolute) max are simple to support and at first glance seem like basic things that we should support, in fact I think they are rather unusual choices (as you say, you would probably only pick min in very particular circumstances / with a lot of prior knowledge about what to expect in that ROI at that time).
supporting a callable seems like a good way to go, but might actually require some refactoring to make a nice API. As you can see, the current implementations are generally lambda flip, data: func(data) (i.e., a function that takes in two params and ignores the first one). If we support a callable it should probably be a callable that takes one argument (the data), not two. This probably means refactoring _pca_flip to be a partialed function, and simplifying the lambda expressions for the other modes.

Aug 25 '23 14:08 drammock

@drammock @agramfort Would it be OK if I start to work on this issue? I am also interested in having this functionality (passing a callable to extract_label_time_course) for my project.

Oct 13 '23 09:10 ctrltz

@drammock @agramfort Would it be OK if I start to work on this issue? I am also interested in having this functionality (passing a callable to extract_label_time_course) for my project.

Go for it! And don't hesitate to ask if you have questions or get stuck on anything.

Oct 13 '23 11:10 drammock

@drammock Thanks! I have looked at the code a bit and wanted to merge all the suggestions with some of my ideas/questions to clarify the details.

All the existing options remain the same, description of max should be fixed in the docs (should I make a separate small PR for this?)
Should I add the center of mass option that was suggested by the author of the issue? I have recently analyzed papers that perform extraction of ROI activity, and averaging, SVD, and centroid/center of mass were the top three options with more or less equal usage. Minimum value indeed sounds like quite an unusual choice.
Regarding supporting mode as a callable: would it make sense to also provide label and src to the custom logic, making it three arguments in total? Having this information, one could apply sign flip if needed by using mne.label_sign_flip(label, src) within the callable or set up weighting schemes that take into account the distance to the center of the ROI.

Nov 18 '23 15:11 ctrltz

All the existing options remain the same, description of max should be fixed in the docs (should I make a separate small PR for this?)

Yes

Should I add the center of mass option

Yes

Regarding supporting mode as a callable ...

I defer to @larsoner on this question (as I'm currently on holiday and can't easily look at the source code, which helps me foresee possible issues)

Nov 21 '23 10:11 drammock

Yes I think it makes sense to pass the label and src to the callable as well

Nov 21 '23 21:11 larsoner

Thanks for the replies @larsoner @drammock! I was trying out possible solutions locally and arrived at several more questions that I would like to discuss.

For the center of mass option, I only later realized that label.center_of_mass has the subjects_dir argument, and in the current implementation I see no way to obtain the suitable value of subjects_dir within mne.source_estimate._gen_extract_label_time_course. Adding subjects_dir as an argument to _gen_extract_label_time_course does not seem nice since other modes like max or mean do not need it at all. If callable was supported, it would be possible to implement the center of mass option as a callable. In that case, it could be added in one of the tutorials/examples but not provided by MNE directly.
I am still trying out different design solutions for the API for the callable. On the one hand, for higher customization it makes sense to provide more information to the callable (e.g., label, src, and maybe even flip). On the other hand, with more arguments the signature of the callable becomes messier especially if some of the arguments are not required (e.g, for taking the minimum value). Trying to combine best of both worlds, I stumbled upon an inspection that determines the number of arguments (e.g., positional ones) of the provided function and thus allows passing a suitable number of input values (see the demo below). Do you think that this is an over-complication or a reasonable way to proceed?

one_arg_fun = lambda data: print(f'Argument(s): {data}')
two_arg_fun = lambda data, label: print(f'Argument(s): {data}, {label}')
three_arg_fun = lambda data, label, src: print(f'Argument(s): {data}, {label}, {src}')
three_arg_fun2 = lambda data, _, src: print(f'Argument(s): {data}, {src}')

def accepts_callable_fun(func: callable):
    args = ['data', 'label', 'src']

    func(*args[:func.__code__.co_argcount])

accepts_callable_fun(one_arg_fun)      # Argument(s): data
accepts_callable_fun(two_arg_fun)      # Argument(s): data, label
accepts_callable_fun(three_arg_fun)    # Argument(s): data, label, src
accepts_callable_fun(three_arg_fun2)   # Argument(s): data, src

I would be very happy to hear your thoughts!

Dec 08 '23 13:12 ctrltz

label.center_of_mass has the subjects_dir argument, and in the current implementation I see no way to obtain the suitable value of subjects_dir within mne.source_estimate._gen_extract_label_time_course.

I think in label.center_of_mass, the subjects_dir is only needed to load the source space surface mesh. In the STC, you already have access to stc.surface(). So some refactoring of _center_of_mass() (the private func, note the leading underscore) to allow it to accept a pre-specified surface (rather than loading it from disk) might be enough to get this working without needing to specify subjects_dir

func(*args[:func.code.co_argcount])