Xfit-style multi-dipole modeling: get timecourses for multiple dipoles

Looking at the git blame for the dipole code, I think I’m addressing this question mostly to @larsoner, but I’d love input from anyone who knows something about for Xfit handles multiple dipoles.

What I’m trying to achieve: Using evoked data, I have created multiple dipoles in MEGIN’s Xfit tool. Now, I would like to take these dipoles and create the corresponding source timecourses for them on single-epoch data.

Problem: Doing this in Xfit would be a cumbersome process, as it doesn’t support epoched data, so I would have to convert a single epoch to an evoked object, load it into Xfit, get the timecourses, convert it back, etc. Not fun. MNE-Python can compute the timecourse given a single dipole position and orientation through the fit_dipole function. However, it cannot do this for multiple dipoles at the same time. It is important to do this at the same time and not dipole-a-dipole, because sensor activity that has been “assigned” to one dipole should not be assigned to the second one. To the best of my knowledge, Xfit achieves this by first whitening the leadfield and evoked data before computing the source timecourses.

Desired solution: I’m trying to modify the fit_dipole function so that multiple values for pos and ori can be supplied corresponding to multiple dipoles. But these functions are long and complicated. Some help would be greatly appreciated.

I wouldn’t do this – multi-dipole fitting in MNE is better handled by the MxNE code. It does nice stuff like split the GOF properly amongst the dipoles used. I would compute a force-fixed forward (normal ori only, not free ori) using a discrete source space created using setup_volume_source_space(pos=dict(rr=rr, nn=nn), ...) where rr and nn are each of shape (2, 3) (rr containing pos and nn containing the ori of your two dipoles of interest), then run mne.inverse_sparse.mixed_norm on it.

1 Like

Thanks for the tip. I’ll give it a try!

It worked! But I had to hack MNE to do it ([TEST] Playing around with fixed dipole orientation in a discrete source space by wmvanvliet · Pull Request #10440 · mne-tools/mne-python · GitHub). Currently, we don’t support volume or discrete forward models with fixed orientation. Is there a strong reason for this?

Here is some example code that currently doesn’t work, but when you apply the PR it does:

import mne
import numpy as np
import matplotlib.pyplot as plt

# Read phantom data
data_path = mne.datasets.brainstorm.bst_phantom_elekta.data_path(verbose=True)
raw = mne.io.read_raw_fif(f'{data_path}/kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif')
events = mne.find_events(raw, 'STI201')

# Extract evoked data for 2 dipoles
epochs = mne.Epochs(raw, events, [1, 17], tmin=-0.1, tmax=0.1, baseline=(-0.1, -0.05))
evoked_dip1 = epochs['1'].average()
evoked_dip2 = epochs['17'].average()

# Create an evoked with both dipole active at the same time
evoked = mne.combine_evoked([evoked_dip1, evoked_dip2], [1, 1])

# Fit two dipoles
cov = mne.compute_covariance(epochs, tmax=-0.05)
bem = mne.make_sphere_model(r0=[0, 0, 0], head_radius=0.08)
dip1, _ = mne.fit_dipole(evoked_dip1.copy().crop(0.035, 0.035), cov, bem)
dip2, _ = mne.fit_dipole(evoked_dip2.copy().crop(0.035, 0.035), cov, bem)

# Collect the dipoles in a single `Dipole` object
dips = mne.Dipole(
    times=np.hstack((dip1.times, dip2.times)),
    pos=np.vstack((dip1.pos, dip2.pos)),
    ori=np.vstack((dip1.ori, dip2.ori)),
    amplitude=np.hstack((dip1.amplitude, dip2.amplitude)),
    gof=np.hstack((dip1.gof, dip2.gof)),
    name=None,
    conf=None,
    khi2=np.hstack((dip1.khi2, dip2.khi2)),
    nfree=dip1.nfree,
)

# Compute dipole timecourses using MxNE
fwd, _ = mne.make_forward_dipole(dips, bem, evoked.info)
stc = mne.inverse_sparse.mixed_norm(evoked, fwd, cov)

# Plot the timecourses
plt.plot(evoked.times, stc.data.T)

Conversation continues on github: [TEST] Playing around with fixed dipole orientation in a discrete source space by wmvanvliet · Pull Request #10440 · mne-tools/mne-python · GitHub