[ENH, WIP] Add Estimation of Significant Connectivity

doc/references.bib

@@ -203,7 +203,7 @@ @article{NolteEtAl2004
   year = {2004}
 }
 
-@INPROCEEDINGS{li_linear_2017,
+@INPROCEEDINGS{LiEtAl2017,
   author = {Li, Adam and Gunnarsdottir, Kristin M. and Inati, Sara and Zaghloul, Kareem and Gale, John and Bulacio, Juan and Martinez-Gonzalez, Jorge and Sarma, Sridevi V.},
   booktitle = {2017 39th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)},
   title = {Linear time-varying model characterizes invasive EEG signals generated from complex epileptic networks},
@@ -227,6 +227,33 @@ @article{ColcloughEtAl2015
   pages = {439--448}
 }
 
+@Article{LiEtAl2021,
+   author = {Adam Li and Chester Huynh and Zachary Fitzgerald and Iahn Cajigas and Damian Brusko and Jonathan Jagid and Angel O. Claudio and Andres M. Kanner and Jennifer Hopp and Stephanie Chen and Jennifer Haagensen and Emily Johnson and William Anderson and Nathan Crone and Sara Inati and Kareem A. Zaghloul and Juan Bulacio and Jorge Gonzalez-Martinez and Sridevi V. Sarma},
+   doi = {10.1038/s41593-021-00901-w},
+   issn = {15461726},
+   issue = {10},
+   journal = {Nature Neuroscience},
+   title = {Neural fragility as an EEG marker of the seizure onset zone},
+   volume = {24},
+   year = {2021},
+}
+
+@article{RecanatesiEtAl2022,
+  title = {Metastable attractors explain the variable timing of stable behavioral action sequences},
+  volume = {110},
+  issn = {08966273},
+  url = {https://linkinghub.elsevier.com/retrieve/pii/S0896627321007790},
+  doi = {10.1016/j.neuron.2021.10.011},
+  language = {en},
+  number = {1},
+  urldate = {2022-03-19},
+  journal = {Neuron},
+  author = {Recanatesi, Stefano and Pereira-Obilinovic, Ulises and Murakami, Masayoshi and Mainen, Zachary and Mazzucato, Luca},
+  month = jan,
+  year = {2022},
+  pages = {139--153.e9},
+}
+
  @article{StamEtAl2012,
   title={Go with the flow: Use of a directed phase lag index (dPLI) to characterize patterns of phase relations in a large-scale model of brain dynamics},
   volume={62},

examples/dynamic/mne_var_connectivity.py

@@ -6,21 +6,21 @@
 ===================================================
 
 Compute a VAR (linear system) model from time-series
-activity :footcite:`li_linear_2017` using a
-continuous iEEG recording.
+activity :footcite:`LiEtAl2017`.
 
 In this example, we will demonstrate how to compute
 a VAR model with different statistical assumptions.
 """
 
 # Authors: Adam Li <adam2392@gmail.com>
+#          Alex Rockhill <aprockhill@mailbox.org>
 #
 # License: BSD (3-clause)
+
 import numpy as np
 
 import mne
-from mne import make_fixed_length_epochs
-from mne_bids import BIDSPath, read_raw_bids
+from mne.datasets import sample
 
 import matplotlib.pyplot as plt
 
@@ -29,80 +29,31 @@
 # %%
 # Load the data
 # -------------
-# Here, we first download an ECoG dataset that was recorded from a patient with
-# epilepsy. To facilitate loading the data, we use `mne-bids
-# <https://mne.tools/mne-bids/>`_.
+# Here, we first download a somatosensory dataset.
 #
 # Then, we will do some basic filtering and preprocessing using MNE-Python.
 
-# paths to mne datasets - sample ECoG
-bids_root = mne.datasets.epilepsy_ecog.data_path()
-
-# first define the BIDS path
-bids_path = BIDSPath(root=bids_root, subject='pt1', session='presurgery',
-                     task='ictal', datatype='ieeg', extension='vhdr')
-
-# Then we'll use it to load in the sample dataset. Here we use a format (iEEG)
-# that is only available in MNE-BIDS 0.7+, so it will emit a warning on
-# versions <= 0.6
-raw = read_raw_bids(bids_path=bids_path, verbose=False)
-
-line_freq = raw.info['line_freq']
-print(f'Data has a power line frequency at {line_freq}.')
-
-# Pick only the ECoG channels, removing the ECG channels
-raw.pick_types(ecog=True)
-
-# Load the data
-raw.load_data()
-
-# Then we remove line frequency interference
-raw.notch_filter(line_freq)
-
-# drop bad channels
-raw.drop_channels(raw.info['bads'])
-
-# %%
-# Crop the data for this example
-# ------------------------------
-#
-# We find the onset time of the seizure and remove all data after that time.
-# In this example, we are only interested in analyzing the interictal
-# (non-seizure) data period.
-#
-# One might be interested in analyzing the seizure period also, which we
-# leave as an exercise for our readers!
-
-# Find the annotated events
-events, event_id = mne.events_from_annotations(raw)
-
-# get sample at which seizure starts
-onset_id = event_id['onset']
-onset_idx = np.argwhere(events[:, 2] == onset_id)
-onset_sample = events[onset_idx, 0].squeeze()
-onset_sec = onset_sample / raw.info['sfreq']
+data_path = sample.data_path()
+raw_fname = data_path / 'MEG' / 'sample' / \
+    'sample_audvis_filt-0-40_raw.fif'
+event_fname = data_path / 'MEG' / 'sample' / \
+    'sample_audvis_filt-0-40_raw-eve.fif'
 
-# remove all data after the seizure onset
-raw = raw.crop(tmin=0, tmax=onset_sec, include_tmax=False)
+# Setup for reading the raw data
+raw = mne.io.read_raw_fif(raw_fname)
+events = mne.read_events(event_fname)
 
-# %%
-# Create Windows of Data (Epochs) Using MNE-Python
-# ------------------------------------------------
-# We have a continuous iEEG snapshot that is about 60 seconds long
-# (after cropping). We would like to estimate a VAR model over a sliding window
-# of 500 milliseconds with a 250 millisecond step size.
-#
-# We can use `mne.make_fixed_length_epochs` to create an Epochs data structure
-# representing this sliding window.
+# Add a bad channel
+raw.info['bads'] += ['MEG 2443']
 
-epochs = make_fixed_length_epochs(raw=raw, duration=0.5, overlap=0.25)
-times = epochs.times
-ch_names = epochs.ch_names
+# Pick MEG gradiometers
+picks = mne.pick_types(raw.info, meg='grad', eeg=False, stim=False, eog=True,
+                       exclude='bads')
 
-print(epochs)
-print(epochs.times)
-print(epochs.event_id)
-print(epochs.events)
+# Create epochs for the visual condition
+event_id, tmin, tmax = 3, -0.2, 1.5  # need a long enough epoch for 5 cycles
+epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+                    baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6))
 
 
 # %%
@@ -114,7 +65,7 @@
 # time-varying linear system.
 
 conn = vector_auto_regression(
-    data=epochs.get_data(), times=times, names=ch_names)
+    data=epochs.get_data(), times=epochs.times, names=epochs.ch_names)
 
 # this returns a connectivity structure over time
 print(conn)
@@ -158,6 +109,49 @@
 im = ax.imshow(rescov, cmap='viridis', aspect='equal', interpolation='none')
 fig.colorbar(im, cax=cax, orientation='horizontal')
 
+# %%
+# Estimate significant connections using a time-shuffled null distribution
+# ------------------------------------------------------------------------
+# We can maintain autocorrelations by shuffling the channels in time relative
+# to one another as explained in :footcite:`RecanatesiEtAl2022`. The pairwise
+# correlations will then be an estimate of connectivity under a null model
+# of uncorrelated neural data.
+
+null_dist = list()
+data = epochs.get_data()
+rng = np.random.default_rng(99)
+for niter in range(10):  # 1000 or more would be reasonable for a real analysis
+    print(f'Computing null connectivity {niter}')
+    for epo_idx in range(data.shape[0]):
+        for ch_idx in range(data.shape[1]):
+            # pick a random starting time for each epoch and channel
+            start_idx = np.round(rng.random() * data.shape[2]).astype(int)
+            data[epo_idx, ch_idx] = np.concatenate(
+                [data[epo_idx, ch_idx, start_idx:],
+                 data[epo_idx, ch_idx, :start_idx]])
+    null_dist.append(vector_auto_regression(
+        data=data, times=epochs.times, names=epochs.ch_names).get_data())
+
+null_dist = np.array(null_dist)
+
+# %%
+# Visualize significant connections over time with animation
+# ----------------------------------------------------------
+# Let's animate over time to visualize the significant connections at each
+# epoch.
+
+con_data = conn.get_data()
+
+# to bonferroni correct across epochs, use the following:
+threshes = np.quantile(abs(null_dist), 1 - (0.05 / con_data.shape[0]),
+                       axis=(0, 1))
+
+# now, plot the connectivity as it changes for each epoch
+n_lines = np.sum(abs(con_data) > threshes, axis=(1, 2))
+fig, ax = plt.subplots(subplot_kw=dict(projection='polar'), figsize=(10, 10))
+anim, ax = conn.plot_circle(n_lines=n_lines, fontsize_names=4,
+                            fig=fig, ax=ax)
+
 # %%
 # Compute one VAR model using all epochs
 # --------------------------------------
@@ -167,7 +161,7 @@
 # epochs.
 
 conn = vector_auto_regression(
-    data=epochs.get_data(), times=times, names=ch_names,
+    data=epochs.get_data(), times=epochs.times, names=epochs.ch_names,
     model='avg-epochs')
 
 # this returns a connectivity structure over time

examples/sensor_connectivity.py

@@ -12,8 +12,6 @@
 #
 # License: BSD (3-clause)
 
-import os.path as op
-
 import mne
 from mne_connectivity import spectral_connectivity_epochs
 from mne.datasets import sample
@@ -24,10 +22,10 @@
 ###############################################################################
 # Set parameters
 data_path = sample.data_path()
-raw_fname = op.join(data_path, 'MEG', 'sample',
-                    'sample_audvis_filt-0-40_raw.fif')
-event_fname = op.join(data_path, 'MEG', 'sample',
-                      'sample_audvis_filt-0-40_raw-eve.fif')
+raw_fname = data_path / 'MEG' / 'sample' / \
+    'sample_audvis_filt-0-40_raw.fif'
+event_fname = data_path / 'MEG' / 'sample' / \
+    'sample_audvis_filt-0-40_raw-eve.fif'
 
 # Setup for reading the raw data
 raw = mne.io.read_raw_fif(raw_fname)

mne_connectivity/base.py

@@ -774,7 +774,7 @@ def rename_nodes(self, mapping):
 
     @copy_function_doc_to_method_doc(plot_connectivity_circle)
     def plot_circle(self, **kwargs):
-        plot_connectivity_circle(
+        return plot_connectivity_circle(
             self.get_data(),
             node_names=self.names,
             indices=self.indices, **kwargs)

mne_connectivity/viz/circle.py

@@ -5,8 +5,10 @@
 #
 # License: Simplified BSD
 
+import numpy as np
 from mne.utils import warn
 from mne.viz.circle import _plot_connectivity_circle
+from matplotlib.animation import FuncAnimation
 
 
 def plot_connectivity_circle(con, node_names, indices=None, n_lines=None,
@@ -105,7 +107,7 @@ def plot_connectivity_circle(con, node_names, indices=None, n_lines=None,
 
     Returns
     -------
-    fig : instance of matplotlib.figure.Figure
+    fig : instance of matplotlib.figure.Figure | instance of matplotlib.animation.FuncAnimation  # noqa E501
         The figure handle.
     ax : instance of matplotlib.projections.polar.PolarAxes
         The subplot handle.
@@ -141,6 +143,40 @@ def plot_connectivity_circle(con, node_names, indices=None, n_lines=None,
                 subplot = (subplot,)
             ax = plt.subplot(*subplot, polar=True)
 
+    if con.ndim == 3:
+        if ax is None:
+            fig, ax = plt.subplots(subplot_kw=dict(projection='polar'),
+                                   facecolor='black')
+        else:
+            fig = ax.figure
+
+        def update_connectivity(i):
+            ax.clear()
+            these_n_lines = n_lines[i] if isinstance(n_lines, np.ndarray) \
+                else n_lines
+            _plot_connectivity_circle(
+                con=con[i], node_names=node_names, indices=indices,
+                n_lines=these_n_lines, node_angles=node_angles,
+                node_width=node_width, node_height=node_height,
+                node_colors=node_colors, facecolor=facecolor,
+                textcolor=textcolor, node_edgecolor=node_edgecolor,
+                linewidth=linewidth, colormap=colormap, vmin=vmin, vmax=vmax,
+                colorbar=(i == 0 and colorbar), title=title,
+                colorbar_size=colorbar_size,
+                colorbar_pos=colorbar_pos, fontsize_title=fontsize_title,
+                fontsize_names=fontsize_names,
+                fontsize_colorbar=fontsize_colorbar, padding=padding, ax=ax,
+                interactive=interactive, node_linewidth=node_linewidth,
+                show=show)
+            # circle is size 10
+            ax.text(3 * np.pi / 4, 25, f'epoch = {i}', color='white',
+                    clip_on=False)
+
+        anim = FuncAnimation(fig, update_connectivity,
+                             frames=con.shape[0],
+                             blit=False, repeat=False)
+        return anim, ax
+
     return _plot_connectivity_circle(
         con=con, node_names=node_names, indices=indices, n_lines=n_lines,
         node_angles=node_angles, node_width=node_width,

mne_connectivity/viz/tests/test_circle.py

@@ -10,6 +10,7 @@
 import matplotlib.pyplot as plt
 from mne.viz import circular_layout
 
+from mne_connectivity import EpochConnectivity
 from mne_connectivity.viz import plot_connectivity_circle
 
 
@@ -92,4 +93,9 @@ def test_plot_connectivity_circle():
                   group_boundaries=[-1])
     pytest.raises(ValueError, circular_layout, label_names, node_order,
                   group_boundaries=[20, 0])
+
+    # test animation
+    conn = EpochConnectivity(data=np.repeat(con.flatten()[None], 10, axis=0),
+                             n_nodes=68)
+    anim, ax = conn.plot_circle(n_lines=3)
     plt.close('all')