Extracting cross-frequency coupling waveforms from rodent LFP data

This example illustrates how to learn univariate atoms on a univariate time-serie. The data is a single LFP channel recorded on a rodent’s striatum 1. Interestingly in this time-serie, the high frequency oscillations around 80 Hz are modulated in amplitude by the low-frequency oscillation around 3 Hz, a phenomenon known as cross-frequency coupling (CFC).

The convolutional sparse coding (CSC) model is able to learn the prototypical waveforms of the signal, on which we can clearly see the CFC.

1

G. Dallérac, M. Graupner, J. Knippenberg, R. C. R. Martinez, T. F. Tavares, L. Tallot, N. El Massioui, A. Verschueren, S. Höhn, J.B. Bertolus, et al. Updating temporal expectancy of an aversive event engages striatal plasticity under amygdala control. Nature Communications, 8:13920, 2017

# Authors: Tom Dupre La Tour <tom.duprelatour@telecom-paristech.fr>
#          Mainak Jas <mainak.jas@telecom-paristech.fr>
#          Umut Simsekli <umut.simsekli@telecom-paristech.fr>
#          Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>
#
# License: BSD (3-clause)

Let us first load the data sample.

import mne
import numpy as np
import matplotlib.pyplot as plt

# sample frequency
sfreq = 350.

# We load the signal. It is an LFP channel recorded on a rodent's striatum.
data = np.load('../rodent_striatum.npy')
print(data.shape)

Out:

(1, 630000)

As the data contains severe artifacts between t=0 and t=100, we use a section not affected by artifacts.

data = data[:, 35000:]

# We also remove the slow drift, which accounts for a lot of variance.
data = mne.filter.filter_data(data, sfreq, 1, None)

# To make the most of parallel computing, we split the data into trials.
data = data.reshape(50, -1)
data /= data.std()

Out:

Setting up high-pass filter at 1 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal highpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Filter length: 1155 samples (3.300 sec)

This sample contains CFC between 3 Hz and 80 Hz. This phenomenon can be described with a comodulogram, computed for instance with the pactools Python library.

from pactools import Comodulogram

comod = Comodulogram(fs=sfreq, low_fq_range=np.arange(0.2, 10.2, 0.2),
                     low_fq_width=2., method='duprelatour')
comod.fit(data)
comod.plot()
plt.show()

Out:

[                                        ] 0% | 0.00 sec | comodulogram: DAR(10, 1)
[                                        ] 2% | 2.85 sec | comodulogram: DAR(10, 1)
[.                                       ] 4% | 4.75 sec | comodulogram: DAR(10, 1)
[..                                      ] 6% | 6.90 sec | comodulogram: DAR(10, 1)
[...                                     ] 8% | 7.79 sec | comodulogram: DAR(10, 1)
[....                                    ] 10% | 8.82 sec | comodulogram: DAR(10, 1)
[....                                    ] 12% | 9.67 sec | comodulogram: DAR(10, 1)
[.....                                   ] 14% | 10.67 sec | comodulogram: DAR(10, 1)
[......                                  ] 16% | 11.51 sec | comodulogram: DAR(10, 1)
[.......                                 ] 18% | 12.38 sec | comodulogram: DAR(10, 1)
[........                                ] 20% | 13.46 sec | comodulogram: DAR(10, 1)
[........                                ] 22% | 14.52 sec | comodulogram: DAR(10, 1)
[.........                               ] 24% | 15.34 sec | comodulogram: DAR(10, 1)
[..........                              ] 26% | 16.11 sec | comodulogram: DAR(10, 1)
[...........                             ] 28% | 16.81 sec | comodulogram: DAR(10, 1)
[............                            ] 30% | 17.57 sec | comodulogram: DAR(10, 1)
[............                            ] 32% | 18.34 sec | comodulogram: DAR(10, 1)
[.............                           ] 34% | 19.11 sec | comodulogram: DAR(10, 1)
[..............                          ] 36% | 19.84 sec | comodulogram: DAR(10, 1)
[...............                         ] 38% | 20.54 sec | comodulogram: DAR(10, 1)
[................                        ] 40% | 21.31 sec | comodulogram: DAR(10, 1)
[................                        ] 42% | 22.08 sec | comodulogram: DAR(10, 1)
[.................                       ] 44% | 22.78 sec | comodulogram: DAR(10, 1)
[..................                      ] 46% | 23.58 sec | comodulogram: DAR(10, 1)
[...................                     ] 48% | 24.35 sec | comodulogram: DAR(10, 1)
[....................                    ] 50% | 25.07 sec | comodulogram: DAR(10, 1)
[....................                    ] 52% | 25.78 sec | comodulogram: DAR(10, 1)
[.....................                   ] 54% | 26.51 sec | comodulogram: DAR(10, 1)
[......................                  ] 56% | 27.27 sec | comodulogram: DAR(10, 1)
[.......................                 ] 58% | 28.00 sec | comodulogram: DAR(10, 1)
[........................                ] 60% | 28.72 sec | comodulogram: DAR(10, 1)
[........................                ] 62% | 29.46 sec | comodulogram: DAR(10, 1)
[.........................               ] 64% | 30.18 sec | comodulogram: DAR(10, 1)
[..........................              ] 66% | 30.91 sec | comodulogram: DAR(10, 1)
[...........................             ] 68% | 31.63 sec | comodulogram: DAR(10, 1)
[............................            ] 70% | 32.35 sec | comodulogram: DAR(10, 1)
[............................            ] 72% | 33.06 sec | comodulogram: DAR(10, 1)
[.............................           ] 74% | 33.77 sec | comodulogram: DAR(10, 1)
[..............................          ] 76% | 34.48 sec | comodulogram: DAR(10, 1)
[...............................         ] 78% | 35.23 sec | comodulogram: DAR(10, 1)
[................................        ] 80% | 35.94 sec | comodulogram: DAR(10, 1)
[................................        ] 82% | 36.65 sec | comodulogram: DAR(10, 1)
[.................................       ] 84% | 37.37 sec | comodulogram: DAR(10, 1)
[..................................      ] 86% | 38.09 sec | comodulogram: DAR(10, 1)
[...................................     ] 88% | 38.80 sec | comodulogram: DAR(10, 1)
[....................................    ] 90% | 39.55 sec | comodulogram: DAR(10, 1)
[....................................    ] 92% | 40.26 sec | comodulogram: DAR(10, 1)
[.....................................   ] 94% | 40.97 sec | comodulogram: DAR(10, 1)
[......................................  ] 96% | 41.68 sec | comodulogram: DAR(10, 1)
[....................................... ] 98% | 42.41 sec | comodulogram: DAR(10, 1)
[........................................] 100% | 43.16 sec | comodulogram: DAR(10, 1)

[........................................] 100% | 43.16 sec | comodulogram: DAR(10, 1)

We fit a CSC model on the data.

from alphacsc import learn_d_z

params = dict(
    n_atoms=3,
    n_times_atom=int(sfreq * 1.0),  # 1000. ms
    reg=5.,
    n_iter=10,
    solver_z='l-bfgs',
    solver_z_kwargs=dict(factr=1e9),
    solver_d_kwargs=dict(factr=1e2),
    random_state=42,
    n_jobs=5,
    verbose=1)

_, _, d_hat, z_hat, _ = learn_d_z(data, **params)

Out:

V_0/10 .........

Plot the temporal patterns. Interestingly, we obtain prototypical waveforms of the signal on which we can clearly see the CFC.

n_atoms, n_times_atom = d_hat.shape
n_columns = min(6, n_atoms)
n_rows = int(np.ceil(n_atoms // n_columns))
figsize = (4 * n_columns, 3 * n_rows)
fig, axes = plt.subplots(n_rows, n_columns, figsize=figsize, sharey=True)
axes = axes.ravel()

for kk in range(n_atoms):
    ax = axes[kk]
    time = np.arange(n_times_atom) / sfreq
    ax.plot(time, d_hat[kk], color='C%d' % kk)
    ax.set_xlim(0, n_times_atom / sfreq)
    ax.set(xlabel='Time (sec)', title="Temporal pattern %d" % kk)
    ax.grid(True)

fig.tight_layout()
plt.show()