Add MNE as option for source recon

[matsvanes]

It would be good if we could use minimum norm estimate (MNE) as an alternative to beamforming. This will especially be useful when we expect highly correlated sources (e.g., in an auditory experiment).

[cgohil8]

A secondary benefit of having this option is we can test osl-dynamics on MNE-style source reconstructed data. I've been getting requests from people that want to apply osl-dynamics to data they've source reconstructed with MNE.

I think what we would add to OSL is basically a wrapper for MNE? With batch processing? Otherwise, surely example scripts/tutorials for source reconstruction are already available in the MNE docs?

[matsvanes]

Agreed. This should be relatively straightforward if we write a wrapper for MNE-Python. See https://mne.tools/stable/auto_tutorials/inverse/30_mne_dspm_loreta.html and https://mne.tools/stable/auto_tutorials/inverse/40_mne_fixed_free.html

[matsvanes]

I think something like this would work @cgohil8: (Note that there is some specific code for covid data involved which is irrelevant. Also the forward model is stored in a different location in newer osl verions)

import mne
import numpy as np
from sys import path
path.append('/ohba/pi/knobre/datasets/covid/scripts/covid-meg/utils')
from directories import dirs, fetch_path

fname = fetch_path('clean', 5, 'mmn')[0]
fname_fwd = os.path.join("/".join(fname.replace('preprocessed', 'source').split("/")[:-2]), fname.split("/")[-1].split("_tsss")[0], "rhino", "coreg", "forward-fwd.fif")

data=mne.io.read_raw(fname, preload=True)
fwd = mne.read_forward_solution(fname_fwd)

data_cov = mne.compute_raw_covariance(data, method="empirical", rank=60)
n_channels = data_cov.data.shape[0]
noise_cov_diag = np.zeros(n_channels)
chantypes=['mag', 'grad']
for type in chantypes:
    # Indices of this channel type
    type_data = data.copy().pick(type, exclude="bads")
    inds = []
    for chan in type_data.info["ch_names"]:
        inds.append(data_cov.ch_names.index(chan))

    # Mean variance of channels of this type
    variance = np.mean(np.diag(data_cov.data)[inds])
    noise_cov_diag[inds] = variance

bads = [b for b in data.info["bads"] if b in data_cov.ch_names]
noise_cov = mne.Covariance(noise_cov_diag, data_cov.ch_names, bads, data.info["projs"], nfree=1e10)

inverse_operator = mne.minimum_norm.make_inverse_operator(
    data.info, fwd, noise_cov)

source = mne.minimum_norm.apply_inverse_raw(data, inverse_operator, method="MNE", lambda2=1)

Note that you could equivalently apply this directly to MNE Evoked data, but you would have to use mne.minimum_norm.apply_inverse.

Also note that there are some options for applying the inverse, including lambda2, which has to be set. If this is the same lambda as in Matlab OSL, this is hardcoded to be 1

The code above returns an MNE VolSourceEstimate, which might be a bit annoying to work with. We can always put the data in a raw array like below

source_info = mne.create_info([f"vertex{i}" for i in range(source.data.shape[0])], data.info['sfreq'])
source = mne.io.RawArray(source.data, source_info)
source._annotations = data.annotations