Brainstorm CTF phantom dataset tutorial

Here we compute the evoked from raw for the Brainstorm CTF phantom tutorial dataset. For comparison, see [1] and:

https://neuroimage.usc.edu/brainstorm/Tutorials/PhantomCtf

References

[1]

François Tadel, Sylvain Baillet, John C. Mosher, Dimitrios Pantazis, and Richard M. Leahy. Brainstorm: a user-friendly application for MEG/EEG analysis. Computational Intelligence and Neuroscience, 2011:1–13, 2011. doi:10.1155/2011/879716.

# Authors: Eric Larson <larson.eric.d@gmail.com>
#
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.

import warnings

import matplotlib.pyplot as plt
import numpy as np

import mne
from mne import fit_dipole
from mne.datasets.brainstorm import bst_phantom_ctf
from mne.io import read_raw_ctf

print(__doc__)

The data were collected with a CTF system at 2400 Hz.

data_path = bst_phantom_ctf.data_path(verbose=True)

# Switch to these to use the higher-SNR data:
# raw_path = op.join(data_path, 'phantom_200uA_20150709_01.ds')
# dip_freq = 7.
raw_path = data_path / "phantom_20uA_20150603_03.ds"
dip_freq = 23.0
erm_path = data_path / "emptyroom_20150709_01.ds"
raw = read_raw_ctf(raw_path, preload=True)

The sinusoidal signal is generated on channel HDAC006, so we can use that to obtain precise timing.

sinusoid, times = raw[raw.ch_names.index("HDAC006-4408")]
plt.figure()
plt.plot(times[times < 1.0], sinusoid.T[times < 1.0])

Let’s create some events using this signal by thresholding the sinusoid.

events = np.where(np.diff(sinusoid > 0.5) > 0)[1] + raw.first_samp
events = np.vstack((events, np.zeros_like(events), np.ones_like(events))).T

The CTF software compensation works reasonably well:

raw.plot()

But here we can get slightly better noise suppression, lower localization bias, and a better dipole goodness of fit with spatio-temporal (tSSS) Maxwell filtering:

raw.apply_gradient_compensation(0)  # must un-do software compensation first
mf_kwargs = dict(origin=(0.0, 0.0, 0.0), st_duration=10.0, st_overlap=True)
raw = mne.preprocessing.maxwell_filter(raw, **mf_kwargs)
raw.plot()

Our choice of tmin and tmax should capture exactly one cycle, so we can make the unusual choice of baselining using the entire epoch when creating our evoked data. We also then crop to a single time point (@t=0) because this is a peak in our signal.

tmin = -0.5 / dip_freq
tmax = -tmin
epochs = mne.Epochs(
    raw, events, event_id=1, tmin=tmin, tmax=tmax, baseline=(None, None)
)
evoked = epochs.average()
evoked.plot(time_unit="s")
evoked.crop(0.0, 0.0)

Let’s use a sphere head geometry model and let’s see the coordinate alignment and the sphere location.

sphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.08)

mne.viz.plot_alignment(
    raw.info, subject="sample", meg="helmet", bem=sphere, dig=True, surfaces=["brain"]
)
del raw, epochs

To do a dipole fit, let’s use the covariance provided by the empty room recording.

raw_erm = read_raw_ctf(erm_path).apply_gradient_compensation(0)
raw_erm = mne.preprocessing.maxwell_filter(raw_erm, coord_frame="meg", **mf_kwargs)
cov = mne.compute_raw_covariance(raw_erm)
del raw_erm

with warnings.catch_warnings(record=True):
    # ignore warning about data rank exceeding that of info (75 > 71)
    warnings.simplefilter("ignore")
    dip, residual = fit_dipole(evoked, cov, sphere, verbose=True)

Compare the actual position with the estimated one.

expected_pos = np.array([18.0, 0.0, 49.0])
diff = np.sqrt(np.sum((dip.pos[0] * 1000 - expected_pos) ** 2))
print(f"Actual pos:     {np.array_str(expected_pos, precision=1)} mm")
print(f"Estimated pos:  {np.array_str(dip.pos[0] * 1000, precision=1)} mm")
print(f"Difference:     {diff:0.1f} mm")
print(f"Amplitude:      {1e9 * dip.amplitude[0]:0.1f} nAm")
print(f"GOF:            {dip.gof[0]:0.1f} %")