{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n=====================================================================\nCompute Phase Slope Index (PSI) in source space for a visual stimulus\n=====================================================================\n\nThis example demonstrates how the Phase Slope Index (PSI) [1] can be computed\nin source space based on single trial dSPM source estimates. In addition,\nthe example shows advanced usage of the connectivity estimation routines\nby first extracting a label time course for each epoch and then combining\nthe label time course with the single trial source estimates to compute the\nconnectivity.\n\nThe result clearly shows how the activity in the visual label precedes more\nwidespread activity (a postivive PSI means the label time course is leading).\n\nReferences\n----------\n[1] Nolte et al. \"Robustly Estimating the Flow Direction of Information in\nComplex Physical Systems\", Physical Review Letters, vol. 100, no. 23,\npp. 1-4, Jun. 2008.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Author: Martin Luessi <mluessi@nmr.mgh.harvard.edu>\n#\n# License: BSD (3-clause)\n\n\nimport numpy as np\n\nimport mne\nfrom mne.datasets import sample\nfrom mne.minimum_norm import read_inverse_operator, apply_inverse_epochs\nfrom mne.connectivity import seed_target_indices, phase_slope_index\n\nprint(__doc__)\n\ndata_path = sample.data_path()\nsubjects_dir = data_path + '/subjects'\nfname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'\nfname_raw = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'\nfname_event = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'\nfname_label = data_path + '/MEG/sample/labels/Vis-lh.label'\n\nevent_id, tmin, tmax = 4, -0.2, 0.3\nmethod = \"dSPM\"  # use dSPM method (could also be MNE or sLORETA)\n\n# Load data\ninverse_operator = read_inverse_operator(fname_inv)\nraw = mne.io.read_raw_fif(fname_raw)\nevents = mne.read_events(fname_event)\n\n# pick MEG channels\npicks = mne.pick_types(raw.info, meg=True, eeg=False, stim=False, eog=True,\n                       exclude='bads')\n\n# Read epochs\nepochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,\n                    baseline=(None, 0), reject=dict(mag=4e-12, grad=4000e-13,\n                                                    eog=150e-6))\n\n# Compute inverse solution and for each epoch. Note that since we are passing\n# the output to both extract_label_time_course and the phase_slope_index\n# functions, we have to use \"return_generator=False\", since it is only possible\n# to iterate over generators once.\nsnr = 1.0  # use lower SNR for single epochs\nlambda2 = 1.0 / snr ** 2\nstcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, method,\n                            pick_ori=\"normal\", return_generator=False)\n\n# Now, we generate seed time series by averaging the activity in the left\n# visual corex\nlabel = mne.read_label(fname_label)\nsrc = inverse_operator['src']  # the source space used\nseed_ts = mne.extract_label_time_course(stcs, label, src, mode='mean_flip')\n\n# Combine the seed time course with the source estimates. There will be a total\n# of 7500 signals:\n# index 0: time course extracted from label\n# index 1..7499: dSPM source space time courses\ncomb_ts = zip(seed_ts, stcs)\n\n# Construct indices to estimate connectivity between the label time course\n# and all source space time courses\nvertices = [src[i]['vertno'] for i in range(2)]\nn_signals_tot = 1 + len(vertices[0]) + len(vertices[1])\n\nindices = seed_target_indices([0], np.arange(1, n_signals_tot))\n\n# Compute the PSI in the frequency range 8Hz..30Hz. We exclude the baseline\n# period from the connectivity estimation\nfmin = 8.\nfmax = 30.\ntmin_con = 0.\nsfreq = raw.info['sfreq']  # the sampling frequency\n\npsi, freqs, times, n_epochs, _ = phase_slope_index(\n    comb_ts, mode='multitaper', indices=indices, sfreq=sfreq,\n    fmin=fmin, fmax=fmax, tmin=tmin_con)\n\n# Generate a SourceEstimate with the PSI. This is simple since we used a single\n# seed (inspect the indices variable to see how the PSI scores are arranged in\n# the output)\npsi_stc = mne.SourceEstimate(psi, vertices=vertices, tmin=0, tstep=1,\n                             subject='sample')\n\n# Now we can visualize the PSI using the plot method. We use a custom colormap\n# to show signed values\nv_max = np.max(np.abs(psi))\nbrain = psi_stc.plot(surface='inflated', hemi='lh',\n                     time_label='Phase Slope Index (PSI)',\n                     subjects_dir=subjects_dir,\n                     clim=dict(kind='percent', pos_lims=(95, 97.5, 100)))\nbrain.show_view('medial')\nbrain.add_label(fname_label, color='green', alpha=0.7)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }
  ], 
  "metadata": {
    "kernelspec": {
      "display_name": "Python 2", 
      "name": "python2", 
      "language": "python"
    }, 
    "language_info": {
      "mimetype": "text/x-python", 
      "nbconvert_exporter": "python", 
      "name": "python", 
      "file_extension": ".py", 
      "version": "2.7.13", 
      "pygments_lexer": "ipython2", 
      "codemirror_mode": {
        "version": 2, 
        "name": "ipython"
      }
    }
  }
}