Allow making inverse solutions from fixed-orientation discrete forward models, enabling multi-dipole modeling

examples/inverse/multi_dipole_model.py

@@ -0,0 +1,156 @@
+# -*- coding: utf-8 -*-
+"""
+=================================================================
+Computing source timecourses with an XFit-like multi-dipole model
+=================================================================
+
+MEGIN's XFit program offers a "guided ECD modeling" interface, where multiple
+dipoles can be fitted by manually selecting subsets of sensors and time ranges.
+Then, source timecourses can be computed using a multi-dipole model. The
+advantage of using a multi-dipole model over fitting each dipole in isolation,
+is that when a source is picked up by more than one dipole, a multi-dipole can
+prevent spatial leakage. This example shows how to build a multi-dipole model
+for estimating source timecourses for evokeds or single epochs.
+
+The XFit program is the recommended approach for obtaining dipole fits because
+it offers a convenient interface for it. These dipoles can then be imported
+into MNE-Python by using the :func:`mne.read_dipole` function for building and
+applying the multi-dipole model. However, this example will also demonstrate
+how to perform guided ECD modeling using only MNE-Python functionality.
+"""
+# Author: Marijn van Vliet <w.m.vanvliet@gmail.com>
+#
+# License: BSD-3-Clause
+
+###############################################################################
+# Importing everything and setting up the data paths for the MNE-Sample
+# dataset.
+import warnings
+
+import mne
+from mne.datasets import sample
+from mne.channels import read_vectorview_selection
+from mne.minimum_norm import (make_inverse_operator, apply_inverse,
+                              apply_inverse_epochs)
+import matplotlib.pyplot as plt
+import numpy as np
+
+data_path = sample.data_path()
+meg_path = data_path / 'MEG' / 'sample'
+raw_fname = meg_path / 'sample_audvis_raw.fif'
+cov_fname = meg_path / 'sample_audvis-shrunk-cov.fif'
+bem_dir = data_path / 'subjects' / 'sample' / 'bem'
+bem_fname =  bem_dir / 'sample-5120-5120-5120-bem-sol.fif'
+
+###############################################################################
+# Read the MEG data from the audvis experiment. Make epochs and evokeds for the
+# left and right auditory conditions.
+raw = mne.io.read_raw_fif(raw_fname)
+raw = raw.pick_types(meg=True, eog=True, stim=True)
+info = raw.info
+
+# Create epochs for auditory events
+events = mne.find_events(raw)
+event_id = dict(right=1, left=2)
+epochs = mne.Epochs(raw, events, event_id,
+                    tmin=-0.1, tmax=0.3, baseline=(None, 0),
+                    reject=dict(mag=4e-12, grad=4000e-13, eog=150e-6))
+
+# Create evokeds for left and right auditory stimulation
+evoked_left = epochs['left'].average()
+evoked_right = epochs['right'].average()
+
+###############################################################################
+# Guided dipole modeling, meaning fitting dipoles to a manually selected subset
+# of sensors as a manually chosen time, can now be performed in MEGINs XFit on
+# the evokeds we computed above. However, it is possible to do it completely
+# in MNE-Python.
+
+# Setup conductor model
+cov = mne.read_cov(cov_fname)
+bem = mne.read_bem_solution(bem_fname)
+
+# Fit two dipoles at t=80ms. The first dipole is fitted using only the sensors
+# on the left side of the helmet. The second dipole is fitted using only the
+# sensors on the right side of the helmet.
+picks_left = read_vectorview_selection('Left', info=info)
+evoked_fit_left = evoked_left.copy().crop(0.08, 0.08)
+evoked_fit_left.pick_channels(picks_left)
+
+picks_right = read_vectorview_selection('Right', info=info)
+evoked_fit_right = evoked_right.copy().crop(0.08, 0.08)
+evoked_fit_right.pick_channels(picks_right)
+
+# Fitting the dipoles with a subset of sensors currently generates warnings
+# regarding the projections. We will ignore these for now.
+with warnings.catch_warnings():
+    warnings.simplefilter("ignore")
+    dip_left, _ = mne.fit_dipole(evoked_fit_left, cov, bem)
+    dip_right, _ = mne.fit_dipole(evoked_fit_right, cov, bem)
+
+###############################################################################
+# Now that we have the location and orientations of the dipoles, compute the
+# full timecourses using MNE, assigning activity to both dipoles at the same
+# time while preventing leakage between the two. We use a very low ``lambda``
+# value to ensure both dipoles are fully used.
+
+fwd, _ = mne.make_forward_dipole([dip_left, dip_right], bem, info)
+
+# Apply MNE inverse
+inv = make_inverse_operator(info, fwd, cov, fixed=True, depth=0)
+stc_left = apply_inverse(evoked_left, inv, method='MNE', lambda2=1E-6)
+stc_right = apply_inverse(evoked_right, inv, method='MNE', lambda2=1E-6)
+
+# Plot the timecourses of the resulting source estimate
+fig, axes = plt.subplots(nrows=2, sharex=True, sharey=True)
+axes[0].plot(stc_left.times, stc_left.data.T)
+axes[0].set_title('Left auditory stimulation')
+axes[0].legend(['Dipole 1', 'Dipole 2'])
+axes[1].plot(stc_right.times, stc_right.data.T)
+axes[1].set_title('Right auditory stimulation')
+axes[1].set_xlabel('Time (s)')
+fig.supylabel('Dipole amplitude')
+
+# We can also fit the timecourses to single epochs. Here, we do it for each
+# experimental condition separately.
+stcs_left = apply_inverse_epochs(epochs['left'], inv, lambda2=1E-6,
+                                 method='MNE')
+stcs_right = apply_inverse_epochs(epochs['right'], inv, lambda2=1E-6,
+                                  method='MNE')
+
+###############################################################################
+# To summarize and visualize the single-epoch dipole amplitudes, we will create
+# a detailed plot of the mean amplitude of the dipoles during different
+# experimental conditions.
+
+# Summarize the single epoch timecourses by computing the mean amplitude from
+# 60-90ms.
+amplitudes_left = []
+amplitudes_right = []
+for stc in stcs_left:
+    amplitudes_left.append(stc.crop(0.06, 0.09).mean().data)
+for stc in stcs_right:
+    amplitudes_right.append(stc.crop(0.06, 0.09).mean().data)
+amplitudes = np.vstack([amplitudes_left, amplitudes_right])
+
+# Visualize the epoch-by-epoch dipole ampltudes in a detailed figure.
+n = len(amplitudes)
+n_left = len(amplitudes_left)
+mean_left = np.mean(amplitudes_left, axis=0)
+mean_right = np.mean(amplitudes_right, axis=0)
+
+fig = plt.figure(figsize=(8, 4))
+plt.scatter(np.arange(n), amplitudes[:, 0], label='Dipole 1')
+plt.scatter(np.arange(n), amplitudes[:, 1], label='Dipole 2')
+transition_point = n_left - 0.5
+plt.plot([0, transition_point], [mean_left[0], mean_left[0]], color='C0')
+plt.plot([0, transition_point], [mean_left[1], mean_left[1]], color='C1')
+plt.plot([transition_point, n], [mean_right[0], mean_right[0]], color='C0')
+plt.plot([transition_point, n], [mean_right[1], mean_right[1]], color='C1')
+plt.axvline(transition_point, color='black')
+plt.xlabel('Epochs')
+plt.ylabel('Dipole amplitude')
+plt.suptitle('Single epoch dipole amplitudes')
+fig.text(0.30, 0.9, 'Left auditory stimulation', ha='center')
+fig.text(0.70, 0.9, 'Right auditory stimulation', ha='center')
+plt.legend()

mne/forward/_make_forward.py

@@ -685,6 +685,10 @@ def make_forward_dipole(dipole, bem, info, trans=None, n_jobs=1, verbose=None):
     -----
     .. versionadded:: 0.12.0
     """
+    if isinstance(dipole, list):
+        from ..dipole import _concatenate_dipoles  # To avoid circular import
+        dipole = _concatenate_dipoles(dipole)
+
     # Make copies to avoid mangling original dipole
     times = dipole.times.copy()
     pos = dipole.pos.copy()

mne/forward/forward.py

@@ -1025,10 +1025,10 @@ def _triage_loose(src, loose, fixed='auto'):
     del fixed
 
     for key, this_loose in loose.items():
-        if key != 'surface' and this_loose != 1:
+        if (key != 'surface' and key != 'discrete') and this_loose != 1:
             raise ValueError(
-                'loose parameter has to be 1 or "auto" for non-surface '
-                f'source spaces, got loose["{key}"] = {this_loose}')
+                'loose parameter has to be 1 or "auto" for non-surface/'
+                f'discrete source spaces, got loose["{key}"] = {this_loose}')
         if not 0 <= this_loose <= 1:
             raise ValueError(
                 f'loose ({key}) must be between 0 and 1, got {this_loose}')

mne/forward/tests/test_make_forward.py

@@ -459,6 +459,12 @@ def test_make_forward_dipole(tmp_path):
     assert isinstance(stc, VolSourceEstimate)
     assert_allclose(stc.times, np.arange(0., 0.003, 0.001))
 
+    # Test passing a list of Dipoles instead of a single Dipole object
+    fwd2, stc2 = make_forward_dipole([dip_even_samp[0], dip_even_samp[1:]],
+                                     sphere, info, trans=fname_trans)
+    assert_array_equal(fwd['sol']['data'], fwd2['sol']['data'])
+    assert_array_equal(stc.data, stc2.data)
+
 
 @testing.requires_testing_data
 def test_make_forward_no_meg(tmp_path):

mne/inverse_sparse/tests/test_gamma_map.py

@@ -130,12 +130,6 @@ def test_gamma_map_vol_sphere():
                                     eeg=False, meg=True)
 
     alpha = 0.5
-    pytest.raises(ValueError, gamma_map, evoked, fwd, cov, alpha,
-                  loose=0, return_residual=False)
-
-    pytest.raises(ValueError, gamma_map, evoked, fwd, cov, alpha,
-                  loose=0.2, return_residual=False)
-
     stc = gamma_map(evoked, fwd, cov, alpha, tol=1e-4,
                     xyz_same_gamma=False, update_mode=2,
                     return_residual=False)

mne/inverse_sparse/tests/test_mxne_inverse.py

@@ -204,16 +204,8 @@ def test_mxne_vol_sphere():
     fwd = mne.make_forward_solution(info, trans=None, src=src,
                                     bem=sphere, eeg=False, meg=True)
 
-    alpha = 80.
-    pytest.raises(ValueError, mixed_norm, evoked, fwd, cov, alpha,
-                  loose=0.0, return_residual=False,
-                  maxit=3, tol=1e-8, active_set_size=10)
-
-    pytest.raises(ValueError, mixed_norm, evoked, fwd, cov, alpha,
-                  loose=0.2, return_residual=False,
-                  maxit=3, tol=1e-8, active_set_size=10)
-
     # irMxNE tests
+    alpha = 80.
     with catch_logging() as log:
         stc = mixed_norm(evoked_l21, fwd, cov, alpha,
                          n_mxne_iter=1, maxit=30, tol=1e-8,

mne/minimum_norm/inverse.py

@@ -1395,6 +1395,23 @@ def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca,
     loose = _triage_loose(forward['src'], loose, fixed)
     del fixed
 
+    # We only support fixed orientations for surface and discrete source
+    # spaces. Not volume or mixed.
+    fixed_inverse = all(v == 0. for v in loose.values())
+    if fixed_inverse:
+        if len(loose) > 1:  # Mixed source space
+            raise ValueError('Computing inverse solutions for mixed source '
+                             'spaces with fixed orientations is not '
+                             'supported.')
+        if 'volume' in loose:
+            raise ValueError('Computing inverse solutions for volume source '
+                             'spaces with fixed orientations is not '
+                             'supported.')
+    if loose.get('volume', 1) < 1:
+        raise ValueError('Computing inverse solutions with restricted '
+                         'orientations (loose < 1) is not supported for '
+                         'volume source spaces.')
+
     # Deal with "depth"
     if exp is not None:
         exp = float(exp)
@@ -1403,10 +1420,10 @@ def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca,
                              f'equal to 0, got {exp}')
         exp = exp or None  # alias 0. -> None
 
-    # put the forward solution in correct orientation
+    # Put the forward solution in correct orientation.
     # (delaying for the case of fixed ori with depth weighting if
     # allow_fixed_depth is True)
-    if loose.get('surface', 1.) == 0. and len(loose) == 1:
+    if fixed_inverse:
         if not is_fixed_orient(forward):
             if allow_fixed_depth:
                 # can convert now
@@ -1424,7 +1441,7 @@ def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca,
                 'Forward operator has fixed orientation and can only '
                 'be used to make a fixed-orientation inverse '
                 'operator.')
-        if loose.get('surface', 1.) < 1. and not forward['surf_ori']:
+        if any(v < 1 for v in loose.values()) and not forward['surf_ori']:
             logger.info('Converting forward solution to surface orientation')
             convert_forward_solution(
                 forward, surf_ori=True, use_cps=use_cps, copy=False)
@@ -1442,7 +1459,7 @@ def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca,
             rank=rank)
 
     # Deal with fixed orientation forward / inverse
-    if loose.get('surface', 1.) == 0. and len(loose) == 1:
+    if fixed_inverse:
         orient_prior = None
         if not is_fixed_orient(forward):
             if depth_prior is not None:
@@ -1454,8 +1471,8 @@ def _prepare_forward(forward, info, noise_cov, fixed, loose, rank, pca,
             convert_forward_solution(
                 forward, surf_ori=True, force_fixed=True,
                 use_cps=use_cps, copy=False)
-    else:
-        if loose.get('surface', 1.) < 1:
+    else:  # Free or loose orientation
+        if any(v < 1 for v in loose.values()):
             assert forward['surf_ori']
         # In theory we could have orient_prior=None for loose=1., but
         # the MNE-C code does not do this

mne/minimum_norm/tests/test_inverse.py

@@ -831,6 +831,25 @@ def test_inverse_operator_volume(evoked, tmp_path):
         apply_inverse(evoked, inv_vol, pick_ori='normal')
 
 
+def test_inverse_operator_discrete(evoked, tmp_path):
+    """Test MNE inverse computation on discrete source space."""
+    # Make discrete source space
+    src = mne.setup_volume_source_space(
+        pos=dict(rr=[[0, 0, 0.1], [0, -0.01, 0.05]],
+                 nn=[[0, 1, 0], [1, 0, 0]]),
+        bem=fname_bem)
+
+    # Perform inverse
+    fwd = mne.make_forward_solution(
+        evoked.info, mne.Transform('head', 'mri'), src, fname_bem)
+    inv = make_inverse_operator(
+        evoked.info, fwd, make_ad_hoc_cov(evoked.info), loose=0, fixed=True,
+        depth=0)
+    stc = apply_inverse(evoked, inv)
+    assert (isinstance(stc, VolSourceEstimate))
+    assert stc.data.shape == (2, len(evoked.times))
+
+
 @pytest.mark.slowtest
 @testing.requires_testing_data
 def test_io_inverse_operator(tmp_path):

mne/utils/docs.py

@@ -801,10 +801,11 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75):
 """
 
 docdict['dipole'] = """
-dipole : instance of Dipole
+dipole : instance of Dipole | list of Dipole
     Dipole object containing position, orientation and amplitude of
     one or more dipoles. Multiple simultaneous dipoles may be defined by
-    assigning them identical times.
+    assigning them identical times. Alternatively, multiple simultaneous
+    dipoles may also be specified as a list of Dipole objects.
 """
 
 docdict['distance'] = """