pactools / pactools

Phase-amplitude coupling (PAC) toolbox

Getting Started with pactools

This package provides tools to estimate phase-amplitude coupling (PAC) in neural time series.

In particular, it implements the driven auto-regressive (DAR) models presented in the reference below [Dupre la Tour et al. 2017].

Read more in the documentation.

Installation

To install pactools, use one of the following two commands:

• Latest stable version:

  pip install pactools

• Development version:

  pip install git+https://github.com/pactools/pactools.git#egg=pactools

To upgrade, use the --upgrade flag provided by pip.

To check if everything worked fine, you can do:

python -c 'import pactools'

and it should not give any error messages.

Phase-amplitude coupling (PAC)

Among the different classes of cross-frequency couplings, phase-amplitude coupling (PAC) - i.e. high frequency activity time-locked to a specific phase of slow frequency oscillations - is by far the most acknowledged. PAC is typically represented with a comodulogram, which shows the strenght of the coupling over a grid of frequencies. Comodulograms can be computed in pactools with more than 10 different methods.

Driven auto-regressive (DAR) models

One of the method is based on driven auto-regressive (DAR) models. As this method models the entire spectrum simultaneously, it avoids the pitfalls related to incorrect filtering or the use of the Hilbert transform on wide-band signals. As the model is probabilistic, it also provides a score of the model goodness of fit via the likelihood, enabling easy and legitimate model selection and parameter comparison; this data-driven feature is unique to such model-based approach.

We recommend using DAR models to estimate PAC in neural time-series. More detail in [Dupre la Tour et al. 2017].

Acknowledgment

This work was supported by the ERC Starting Grant SLAB ERC-YStG-676943 to Alexandre Gramfort, the ERC Starting Grant MindTime ERC-YStG-263584 to Virginie van Wassenhove, the ANR-16-CE37-0004-04 AutoTime to Valerie Doyere and Virginie van Wassenhove, and the Paris-Saclay IDEX NoTime to Valerie Doyere, Alexandre Gramfort and Virginie van Wassenhove,

Cite this work

If you use this code in your project, please cite [Dupre la Tour et al. 2017]:

@article{duprelatour2017nonlinear,
    author = {Dupr{\'e} la Tour, Tom and Tallot, Lucille and Grabot, Laetitia and Doy{\`e}re, Val{\'e}rie and van Wassenhove, Virginie and Grenier, Yves and Gramfort, Alexandre},
    journal = {PLOS Computational Biology},
    publisher = {Public Library of Science},
    title = {Non-linear auto-regressive models for cross-frequency coupling in neural time series},
    year = {2017},
    month = {12},
    volume = {13},
    url = {https://doi.org/10.1371/journal.pcbi.1005893},
    pages = {1-32},
    number = {12},
    doi = {10.1371/journal.pcbi.1005893}
}

pactools's Issues

Deprecation warning

@TomDLT wonderful library ;)

I have an issue with warnings. When I run the code in Python 3.6 I get a deprecation warning with n_jobs which is polluting my output. Any idea how I can get rid of it?

/autofs/space/meghnn_001/users/mjas/anaconda3/envs/mne/lib/python3.6/site-packages/_pytest/assertion/rewrite.py:8: DeprecationWarning: the imp module is deprecated in favour of importlib; see the module's documentation for alternative uses

Code to reproduce:

import numpy as np
import matplotlib.pyplot as plt

from pactools import Comodulogram
from pactools import simulate_pac
fs = 200.  # Hz
high_fq = 50.0  # Hz
low_fq = 5.0  # Hz
low_fq_width = 1.0  # Hz

n_points = 1000
noise_level = 0.4

signal = simulate_pac(n_points=n_points, fs=fs, high_fq=high_fq, low_fq=low_fq,
                      low_fq_width=low_fq_width, noise_level=noise_level,
                      random_state=0)
low_fq_range = np.linspace(1, 10, 50)
method = 'duprelatour'  # or 'tort', 'ozkurt', 'penny', 'colgin', ...
n_surrogates = 200
n_jobs = 2
estimator = Comodulogram(fs=fs, low_fq_range=low_fq_range,
                         low_fq_width=low_fq_width, method=method,
                         n_surrogates=n_surrogates, progress_bar=True,
                         n_jobs=n_jobs)
estimator.fit(signal)

The higher the n_jobs, the more I see the warnings. Tried to suppress it using warnings.simplefilter but it doesn't seem to work:

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    from pactools import Comodulogram

Band pass in delay estimator

If I understand the delay estimator correctly it estimates the delay on the entire high frequency part of the PAC (dar model) (assuming the driver is low freq and the pac high frequencies e.g. theta gamma). Would it be possible/easy to implement a way to define which high frequency band to estimate delay for?

I'm asking because I typically get PAC looking like: (from entorhinal cortex L2)

When estimating delay, I'm not sure if it corresponds to the high frequency blob or the low frequency blob.

[Request] Save comodulogram

I'm currently working on a project and I plan to use this tool to compute PAC on several subjects. Since it takes quite some time, I would like to save the comodulograms on files, to be able to do the analysis later.

Are you interested in such a feature? I will do it in any case, but if you are interested, we can discuss a proper way and add it to your tool.

[BUG] Save estimator hdf5 with method as DAR

See minimally reproducible example:

import numpy as np
import matplotlib.pyplot as plt

from pactools import Comodulogram
from pactools import simulate_pac
from pactools.dar_model import DAR, extract_driver
fs = 200.  # Hz
high_fq = 50.0  # Hz
low_fq = 5.0  # Hz
low_fq_width = 1.0  # Hz

n_points = 10000
noise_level = 0.4

signal = simulate_pac(n_points=n_points, fs=fs, high_fq=high_fq, low_fq=low_fq,
                      low_fq_width=low_fq_width, noise_level=noise_level,
                      random_state=0)
dar = DAR(ordar=20, ordriv=2, criterion='bic')
estimator = Comodulogram(fs=fs, low_fq_range=low_fq_range,
                             low_fq_width=low_fq_width, method=dar,
                             progress_bar=False)
estimator.fit(signal)
estimator.save('test.hdf5')

Gives

/Users/alexrockhill/software/anaconda3/envs/swannlab/lib/python3.7/importlib/_bootstrap.py:219: RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 192 from C header, got 216 from PyObject
  return f(*args, **kwds)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/alexrockhill/software/pactools/pactools/comodulogram.py", line 591, in save
    slash='replace')
  File "/Users/alexrockhill/software/anaconda3/envs/swannlab/lib/python3.7/site-packages/mne-0.21.dev0-py3.7.egg/mne/externals/h5io/_h5io.py", line 115, in write_hdf5
    use_json=use_json)
  File "/Users/alexrockhill/software/anaconda3/envs/swannlab/lib/python3.7/site-packages/mne-0.21.dev0-py3.7.egg/mne/externals/h5io/_h5io.py", line 151, in _triage_write
    where + '["%s"]' % key, cleanup_data=cleanup_data, slash=slash)
  File "/Users/alexrockhill/software/anaconda3/envs/swannlab/lib/python3.7/site-packages/mne-0.21.dev0-py3.7.egg/mne/externals/h5io/_h5io.py", line 232, in _triage_write
    raise TypeError(err_str)
TypeError: unsupported type <class 'pactools.dar_model.dar.DAR'> (in <class 'dict'>["method"]) 

saving Comodulogram with BIC selection

When estimating a dar with BIC selection to get optimal hyperparameters, giving the dar as a method in Comodulogram I get the following error

TypeError: unsupported type <class 'pactools.dar_model.dar.DAR'> (in <class 'dict'>["method"])

singular matrices

Hi!

Thank you so much for making this toolset! Finding the method you seek already implemented, well written and well documented is a very joyful experience!

I get LinalgError on some of my datasets, because of singular matrices, have you considered adding a way to estimate the inverse of singular matrices, like a least squares estimation?

Calculation of PAC with TORT method

I am having trouble using the tort method to calculate the PAC and generate a comodulogram for the openneuro datasets (https://openneuro.org/datasets/ds002778/versions/1.0.3). I am calculating the PAC for phase frequency of 13-30hz , bandwidth of 2 and amplitude frequency of 50-100hz bandwidth of 4 for C3 and C4 channel extracted from the EEG dataset. The average PAC values for C3 and C4 together that I get are very small compared to what has already been published using this dataset. For example, file number pd10 , there should be a high PAC value (avg PAC for c3 and c4) for the 'off' session around 0.001 but I am getting around 2.06 e-06. For file number pd14, there should be a really low PAC of 0.00001 but I am getting an even lower value of 1.37 e-06. Also,the comodulogram that is generated is fuzzy even when I use a higher resolution parameter. I am not sure where the issue is.

Here is the code I have used : https://gist.github.com/apoorva6262/6007f348976347272cfb41e670c3e57d
I can share the preprocessed files with you which might be quicker to run. Please let me know.

Additionally, I was going through the source code for comodulogram , and I was not sure if my signal (data1[0]) when I am calling to fit , would be’ high_sig’ or ‘low_sig’. Or is that understood based on the phase and amp frequency ranges that I have listed ?

For the Kullback-Leibler divergence of the distribution , the published paper used :
MI=(log(nbin)-(-sum((MeanAmp/sum(MeanAmp)).*log((MeanAmp/sum(MeanAmp))))))/log(nbin);

But in your source code it shows this , wondering if that would make a difference?
Normalizedmeanamp=MeanAmp/sum(MeanAmp);
Entropy= -sum(Normalizedmeanamp.log(Normalizedmeanamplog(nbin)));
MI= (log(nbin)-(Entropy))/log(nbin);

Any help would be appreciated. Thank you so much.

[BUG] Tort `amplitude_dist` index error

Here is a minimally reproducible example of for some reason it looks like the tort bins don't align.

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/alexrockhill/software/pactools/pactools/comodulogram.py", line 252, in fit
    this_mask, filtered_low_2)
  File "/Users/alexrockhill/software/pactools/pactools/comodulogram.py", line 663, in _comodulogram
    ax_special=estimator.ax_special) for j in range(n_high))
  File "/Users/alexrockhill/software/anaconda3/envs/swannlab/lib/python3.7/site-packages/joblib/parallel.py", line 1004, in __call__
    if self.dispatch_one_batch(iterator):
  File "/Users/alexrockhill/software/anaconda3/envs/swannlab/lib/python3.7/site-packages/joblib/parallel.py", line 835, in dispatch_one_batch
    self._dispatch(tasks)
  File "/Users/alexrockhill/software/anaconda3/envs/swannlab/lib/python3.7/site-packages/joblib/parallel.py", line 754, in _dispatch
    job = self._backend.apply_async(batch, callback=cb)
  File "/Users/alexrockhill/software/anaconda3/envs/swannlab/lib/python3.7/site-packages/joblib/_parallel_backends.py", line 209, in apply_async
    result = ImmediateResult(func)
  File "/Users/alexrockhill/software/anaconda3/envs/swannlab/lib/python3.7/site-packages/joblib/_parallel_backends.py", line 590, in __init__
    self.results = batch()
  File "/Users/alexrockhill/software/anaconda3/envs/swannlab/lib/python3.7/site-packages/joblib/parallel.py", line 256, in __call__
    for func, args, kwargs in self.items]
  File "/Users/alexrockhill/software/anaconda3/envs/swannlab/lib/python3.7/site-packages/joblib/parallel.py", line 256, in <listcomp>
    for func, args, kwargs in self.items]
  File "/Users/alexrockhill/software/pactools/pactools/comodulogram.py", line 676, in _loop_over_shifts
    return [func(shift=sh, **kwargs) for sh in shifts]
  File "/Users/alexrockhill/software/pactools/pactools/comodulogram.py", line 676, in <listcomp>
    return [func(shift=sh, **kwargs) for sh in shifts]
  File "/Users/alexrockhill/software/pactools/pactools/comodulogram.py", line 725, in _one_modulation_index
    amplitude_dist[b] = np.mean(selection)
IndexError: index 18 is out of bounds for axis 0 with size 18

import mne
import numpy as np
import os.path as op
import matplotlib.pyplot as plt

from pactools import simulate_pac, raw_to_mask, Comodulogram, MaskIterator


n_events = 100
mu = 1.  # mean onset of PAC in seconds
sigma = 0.25  # standard deviation of onset of PAC in seconds
trial_len = 2.  # len of the simulated trial in seconds
first_samp = 5  # seconds before the first sample and after the last

fs = 200.  # Hz
high_fq = 50.0  # Hz
low_fq = 3.0  # Hz
low_fq_width = 2.0  # Hz

n_points = int(trial_len * fs)
noise_level = 0.4


def gaussian1d(array, mu, sigma):
    return np.exp(-0.5 * ((array - mu) / sigma) ** 2)


signal = np.zeros((2, int(n_points * n_events + 2 * first_samp * fs)))
events = np.zeros((n_events, 3))
events[:, 0] = np.arange((first_samp + mu) * fs,
                         n_points * n_events + (first_samp + mu) * fs,
                         trial_len * fs, dtype=int)
events[:, 2] = np.ones((n_events))
mod_fun = gaussian1d(np.arange(n_points), mu * fs, sigma * fs)
for i in range(n_events):
    signal_no_pac = simulate_pac(n_points=n_points, fs=fs,
                                 high_fq=high_fq, low_fq=low_fq,
                                 low_fq_width=low_fq_width,
                                 noise_level=1.0, random_state=i)
    signal_pac = simulate_pac(n_points=n_points, fs=fs,
                              high_fq=high_fq, low_fq=low_fq,
                              low_fq_width=low_fq_width,
                              noise_level=noise_level, random_state=i)
    signal[0, i * n_points:(i + 1) * n_points] = \
        signal_pac * mod_fun + signal_no_pac * (1 - mod_fun)

info = mne.create_info(['Ch1', 'STI 014'], fs, ['eeg', 'stim'])
raw = mne.io.RawArray(signal, info)

tmin, tmax = mu - 3 * sigma, mu + 3 * sigma
ixs = (0, 0)
low_sig, high_sig, mask = raw_to_mask(raw, ixs=ixs, events=events, tmin=tmin,
                                      tmax=tmax)
estimator = Comodulogram(fs=raw.info['sfreq'],
                         low_fq_range=np.linspace(1, 10, 20), low_fq_width=2.,
                         method='tort', progress_bar=True)
# compute the comodulogram
estimator.fit(low_sig, high_sig, mask)
# plot the results
estimator.plot(tight_layout=False)
plt.show()

[ENH] Add a decimate keyword argument to Comodulogram

It would be nice if the Comodulogram instance could decimate the signal without causing jitter in the events so that a dataset that was sampled at a high frequency could be computed at a more manageable sampling frequency.

I'm not sure that this is possible but maybe someone else has an idea.

Masking clarification

Thanks for a great set of PAC tools. I'm currently comparing the results with my MATLAB code. I have a couple of questions from the plot_mne_epoch_masking tutorial:

• Do you mask the data after band-pass filtering?
• Do you calculate PAC on single trials, then average... or use the continuously masked data?

(Question) Interpreting phase

Hello,

Could you clarify how to interpret the driver phase? In the driver waveform, I believe 0 degrees corresponds to the peak, 180d is the trough, -90d is the rising zero crossing, and 90d is the falling zero crossing -- is this correct?

If so, why is this the case, instead of 90d corresponding to the peak, and -90d corresponding to the trough, as phase is normally interpreted for a sine wave? I first thought this was due to the Hilbert transform, but I now understand the default bandpass filter does not use the Hilbert transform.

Thank you!

Using pactools for resting state EEG

Pardon me,this is more a research question than usage. I am trying to find differences in cfc between healthy and alzheimer's with little success.I am curious to know if phase amplitude coupling is useful for resting state EEG (no events at all).Any usecases of pactools being used on rest eeg ?

Another important question i got is: Is it better to apply pac on source level (after LORETA source estimation) or is it applied at the sensor (eeg electrode) level too ?

hard to reach raw results

I like the plotting methods, but sometimes you want to reach the result directly. Ended up poking around in the source code to figure out which variable contained the comodulogram. Especially, how to get the spec from the dar is not very clear.

[General usage question-Feature Engineering]How to obtain unique PAC values for a group of patients in the 1020 electrode system?

Hi
Noobie question here - I have a large number of patients whose EEG is recorded using 19 electrodes (single trial, eyes closed, resting). I would like to know, if there is a way to obtain the PAC values for a patient (say, by each electrode) ? I am trying to incorporate CFC PAC as a possible feature in my EEG dataset (after feature engineering/building from raw eeg data using various non linear and linear techniques).
The comod_ attribute in the Comodulogram() object provides a (low_freq_range, high_freq_range) 2d array for each electrode for each patient, as per my present code/understanding. I'd like to know a way to consolidate the comodulogram values so as to represent that with a single feature (per electrode) ? Can you please provide an example of how this can be done ?

Question about default filter design

Hi, I just had a few questions about the default filter design (more of signal processing questions than directly about the code).

1. How does bandwidth translate to the filter? From the code, it's clear the input bandwidth effects the order of the filter, with larger bandwidths corresponding to higher order (and therefore wider) filters. Was this done with any specific benchmark in mind (something like, the order is set so that the stopband of the filter is equal to the center frequency + / - the bandwidth)?

2. What was the reasoning for choosing a filter based off a blackman window? More specifically, here's the filter response of a filter from pactools with a center frequency of 50Hz and a bandwidth of 20Hz:
   

The filter has a gradual slope, so it looks like it still includes a lot of frequencies above or below 40-60Hz. Why not use a filter with a sharper transition?

3. On a related note, do you have any criteria for excluding certain frequency pairs where there is overlap between the low and high frequencies? For the filter above, my first instinct would be to exclude any comparisons between that 50Hz signal and anything above around 16Hz, where the filter response becomes bumpy -- but maybe at frequencies where the response is say -30 dB or less, the frequency can reasonably be considered filtered out.

Thanks so much!

Tort Phase Bins

Hi,

Thanks for a great toolbox!
I wanted to clarify over what phase range the MI is calculated over for the Tort method? It looks like the phase bins range from -pi to pi (phase_bins = np.linspace(-np.pi, np.pi, n_bins)), but when plotting the Phase-Amplitude, the phase range for the x-axis is set to be phase_bins = 0.5 * (phase_bins[:-1] + phase_bins[1:]) / np.pi * 180 (-180 to 180, or -pi / 2 to pi / 2).

[ENH] Use simulated data instead of MNE sample data for I/O demo

Since I already wrote most of how to do that here: #22, I could change a few lines around so that you make an MNE object with PAC, save it to a temp directory, load it in and compute PAC. IMO that would make a better example.

This is in reference to this demo: https://pactools.github.io/auto_examples/plot_mne_epoch_masking.html.

Implement cross-bispectrum

Currently, bispectrum-based metrics can only be calculated on a single channel.
Adding a cross-channel implementation would be valuable.

[ENH] Add an event-related PAC simulation

Maybe @TomDLT would be interested in reviewing, if I were to work on a PR to simulate PAC in an event-related manner with a block design if that is something of interest. Similar to the 2019 JNeurosci paper but I was expecting instead of an a priori region of interest (0.4 to 1.2 seconds between two button presses for instance) a search/method to determine the region of interest. Or maybe, even more ambitiously, a measure of instantaneous PAC using many similar trials using something possibly like the Morlet approach if that would have any chance of working.

Where did randomisation procedure happen and what the purpose of it?

First of all, this is a great work and I found it works very well with simulated data. I would like to explore more with my real-life data. However, everytime I calculated the PAC with this toolbox with method='duprelatour', I got different result. I understand that surrogate may have happened but since I would like to control this random process therefore I would rather shift the phase by my own function instead of using randomisation procedure.
However, even if I use the parameter like this:

comod = Comodulogram(fs=epochs.info['sfreq'], low_fq_range=fq_range[0], high_fq_range=fq_range[1], 
                     n_surrogates=0, n_jobs=None,
                     progress_bar=True, 
                     minimum_shift=1,
                     method='duprelatour')

I still got different results every time. I am wondering of which function did the randomisation happen?
Thanks.

binary incompatibility warning

Dear Developers,

I have installed the pactools using pip command. My python version is 3.8.5 under win 10, and conda version 4.10.3.
After importing the library I got the following warnings: