Note

Click here to download the full example code

CSC to learn LFP spiking atoms#

Here, we show how CSC can be used to learn spiking atoms from Local Field Potential (LFP) data [1].

[1] Hitziger, Sebastian, et al.

Adaptive Waveform Learning: A Framework for Modeling Variability in Neurophysiological Signals. IEEE Transactions on Signal Processing (2017).

First, let us fetch the data (~14 MB)

import os
from mne.datasets import fetch_dataset
from mne.utils import get_config

url = ('https://github.com/hitziger/AWL/raw/master/Experiments/data/'
       'LFP_data_contiguous_1250_Hz.mat')

folder_name = "LFP"
archive_name = "LFP_data_contiguous_1250_Hz.mat"

fname = fetch_dataset(
    {"dataset_name": "LFP_data",
     "url": url,
     "archive_name": archive_name,
     "folder_name": folder_name,
     "hash": None
     },
    path=None,
    force_update=False
)

fname = os.path.join(fname, archive_name)

Out:

Using default location ~/mne_data for LFP_data...
Creating ~/mne_data
Downloading file 'LFP_data_contiguous_1250_Hz.mat' from 'https://github.com/hitziger/AWL/raw/master/Experiments/data/LFP_data_contiguous_1250_Hz.mat' to '/github/home/mne_data/LFP'.

  0%|                                              | 0.00/14.4M [00:00<?, ?B/s]
 46%|█████████████████                    | 6.64M/14.4M [00:00<00:00, 66.4MB/s]
 93%|██████████████████████████████████▎  | 13.3M/14.4M [00:00<00:00, 66.8MB/s]
  0%|                                              | 0.00/14.4M [00:00<?, ?B/s]
100%|█████████████████████████████████████| 14.4M/14.4M [00:00<00:00, 28.7GB/s]

It is a mat file, so we use scipy to load it

from scipy import io

data = io.loadmat(fname)
X, sfreq = data['X'].T, float(data['sfreq'])

And now let us look at the data

import numpy as np
import matplotlib.pyplot as plt

start, stop = 11000, 15000
times = np.arange(start, stop) / sfreq
plt.plot(times, X[0, start:stop], color='b')
plt.xlabel('Time (s)')
plt.ylabel(r'$\mu$ V')
plt.xlim([9., 12.])
plot lfp data

Out:

(9.0, 12.0)

and filter it using a convenient function from MNE. This will remove low frequency drifts, but we keep the high frequencies

from mne.filter import filter_data
X = filter_data(
    X.astype(np.float64), sfreq, l_freq=1, h_freq=None, fir_design='firwin')

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: 4125 samples (3.300 sec)

Now, we define the parameters of our model.

reg = 6.0
n_times = 2500
n_times_atom = 350
n_trials = 100
n_atoms = 3
n_iter = 60

Let’s stick to one random state for now, but if you want to learn how to select the random state, consult this example.

random_state = 10

Now, we epoch the trials

overlap = 0
starts = np.arange(0, X.shape[1] - n_times, n_times - overlap)
stops = np.arange(n_times, X.shape[1], n_times - overlap)

X_new = []
for idx, (start, stop) in enumerate(zip(starts, stops)):
    if idx >= n_trials:
        break
    X_new.append(X[0, start:stop])
X_new = np.vstack(X_new)
del X

We remove the mean and scale to unit variance.

X_new -= np.mean(X_new)
X_new /= np.std(X_new)

The convolutions can result in edge artifacts at the edges of the trials. Therefore, we discount the contributions from the edges by windowing the trials.

from numpy import hamming
X_new *= hamming(n_times)[None, :]

Of course, in a data-limited setting we want to use as much of the data as possible. If this is the case, you can set overlap to non-zero (for example half the epoch length).

Now, we run regular CSC since the trials are not too noisy

from alphacsc import learn_d_z
pobj, times, d_hat, z_hat, reg = learn_d_z(X_new, n_atoms, n_times_atom,
                                           reg=reg, n_iter=n_iter,
                                           random_state=random_state, n_jobs=1)

Out:

Coordinate descent loop 0 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 51134.32096508
Using method dual for projection
[seed 10] Objective (d) 45299.80903125
Coordinate descent loop 1 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 26399.47853871
Using method dual for projection
[seed 10] Objective (d) 25948.80915409
Coordinate descent loop 2 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 23777.89064427
Using method dual for projection
[seed 10] Objective (d) 23704.90593344
Coordinate descent loop 3 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 22233.60227587
Using method dual for projection
[seed 10] Objective (d) 22221.01870097
Coordinate descent loop 4 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 21202.43876006
Using method dual for projection
[seed 10] Objective (d) 21169.49202496
Coordinate descent loop 5 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 20579.63611160
Using method dual for projection
[seed 10] Objective (d) 20618.98107838
Coordinate descent loop 6 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 20283.55110757
Using method dual for projection
[seed 10] Objective (d) 20210.75292119
Coordinate descent loop 7 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 19971.53116238
Using method dual for projection
[seed 10] Objective (d) 19970.52325693
Coordinate descent loop 8 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 19708.44554814
Using method dual for projection
[seed 10] Objective (d) 19691.20196937
Coordinate descent loop 9 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 19462.43966049
Using method dual for projection
[seed 10] Objective (d) 19467.01239142
Coordinate descent loop 10 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 19233.45092821
Using method dual for projection
[seed 10] Objective (d) 19214.12766228
Coordinate descent loop 11 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 19085.64469638
Using method dual for projection
[seed 10] Objective (d) 19084.13170902
Coordinate descent loop 12 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18945.61954180
Using method dual for projection
[seed 10] Objective (d) 18980.31578130
Coordinate descent loop 13 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18876.43985522
Using method dual for projection
[seed 10] Objective (d) 18836.38353365
Coordinate descent loop 14 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18720.74584344
Using method dual for projection
[seed 10] Objective (d) 18721.54845149
Coordinate descent loop 15 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18623.64002131
Using method dual for projection
[seed 10] Objective (d) 18621.18904484
Coordinate descent loop 16 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18547.81611677
Using method dual for projection
[seed 10] Objective (d) 18545.24281149
Coordinate descent loop 17 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18471.88752076
Using method dual for projection
[seed 10] Objective (d) 18470.42055102
Coordinate descent loop 18 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18393.84202336
Using method dual for projection
[seed 10] Objective (d) 18391.78226334
Coordinate descent loop 19 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18331.73169086
Using method dual for projection
[seed 10] Objective (d) 18328.84451729
Coordinate descent loop 20 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18274.15271778
Using method dual for projection
[seed 10] Objective (d) 18297.40702007
Coordinate descent loop 21 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18257.44912760
Using method dual for projection
[seed 10] Objective (d) 18276.17808963
Coordinate descent loop 22 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18243.49274421
Using method dual for projection
[seed 10] Objective (d) 18261.00618790
Coordinate descent loop 23 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18234.36110400
Using method dual for projection
[seed 10] Objective (d) 18208.94505630
Coordinate descent loop 24 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18179.04646062
Using method dual for projection
[seed 10] Objective (d) 18179.60670430
Coordinate descent loop 25 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18156.84098639
Using method dual for projection
[seed 10] Objective (d) 18166.58627562
Coordinate descent loop 26 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18145.85002356
Using method dual for projection
[seed 10] Objective (d) 18152.93057120
Coordinate descent loop 27 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18131.29820966
Using method dual for projection
[seed 10] Objective (d) 18131.16931585
Coordinate descent loop 28 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18119.63495602
Using method dual for projection
[seed 10] Objective (d) 18111.53150159
Coordinate descent loop 29 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18095.94336614
Using method dual for projection
[seed 10] Objective (d) 18106.50704724
Coordinate descent loop 30 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18098.05237401
Using method dual for projection
[seed 10] Objective (d) 18097.82987219
Coordinate descent loop 31 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18084.15728858
Using method dual for projection
[seed 10] Objective (d) 18025.10056909
Coordinate descent loop 32 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18015.48732444
Using method dual for projection
[seed 10] Objective (d) 18015.15575745
Coordinate descent loop 33 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18009.09807550
Using method dual for projection
[seed 10] Objective (d) 18010.63887828
Coordinate descent loop 34 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18004.88799181
Using method dual for projection
[seed 10] Objective (d) 18058.43295042
Coordinate descent loop 35 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18049.18814404
Using method dual for projection
[seed 10] Objective (d) 17997.43386062
Coordinate descent loop 36 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17990.54845736
Using method dual for projection
[seed 10] Objective (d) 17990.87720230
Coordinate descent loop 37 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17985.96104110
Using method dual for projection
[seed 10] Objective (d) 17988.72476199
Coordinate descent loop 38 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17984.75059699
Using method dual for projection
[seed 10] Objective (d) 18049.42022273
Coordinate descent loop 39 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 18046.54425058
Using method dual for projection
[seed 10] Objective (d) 17984.35856862
Coordinate descent loop 40 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17980.96880430
Using method dual for projection
[seed 10] Objective (d) 17980.80669037
Coordinate descent loop 41 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17972.69296219
Using method dual for projection
[seed 10] Objective (d) 17972.04008927
Coordinate descent loop 42 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17966.98480639
Using method dual for projection
[seed 10] Objective (d) 17966.82982420
Coordinate descent loop 43 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17962.43826374
Using method dual for projection
[seed 10] Objective (d) 17977.53550322
Coordinate descent loop 44 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17975.03587533
Using method dual for projection
[seed 10] Objective (d) 17966.16026575
Coordinate descent loop 45 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17963.76015273
Using method dual for projection
[seed 10] Objective (d) 17978.91055034
Coordinate descent loop 46 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17976.75050954
Using method dual for projection
[seed 10] Objective (d) 17972.49942963
Coordinate descent loop 47 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17970.48291179
Using method dual for projection
[seed 10] Objective (d) 17973.23422992
Coordinate descent loop 48 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17970.92478592
Using method dual for projection
[seed 10] Objective (d) 17963.82948038
Coordinate descent loop 49 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17961.84373701
Using method dual for projection
[seed 10] Objective (d) 17971.86050776
Coordinate descent loop 50 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17967.92414417
Using method dual for projection
[seed 10] Objective (d) 17974.11292268
Coordinate descent loop 51 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17971.19344934
Using method dual for projection
[seed 10] Objective (d) 17976.59916952
Coordinate descent loop 52 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17974.68790064
Using method dual for projection
[seed 10] Objective (d) 17972.97115106
Coordinate descent loop 53 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17971.74333296
Using method dual for projection
[seed 10] Objective (d) 17971.21362948
Coordinate descent loop 54 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17970.14057303
Using method dual for projection
[seed 10] Objective (d) 17963.95975857
Coordinate descent loop 55 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17962.82815869
Using method dual for projection
[seed 10] Objective (d) 17964.32024402
Coordinate descent loop 56 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17963.16483309
Using method dual for projection
[seed 10] Objective (d) 17960.62470416
Coordinate descent loop 57 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17959.04029128
Using method dual for projection
[seed 10] Objective (d) 17956.00852062
Coordinate descent loop 58 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17954.27151859
Using method dual for projection
[seed 10] Objective (d) 17956.37784603
Coordinate descent loop 59 / 60 [n_jobs=1]
[seed 10] Objective (z_hat) : 17955.30068453
Using method dual for projection
[seed 10] Objective (d) 17945.26047536

Let’s look at the atoms now.

plt.figure()
plt.plot(d_hat.T)
plt.show()
plot lfp data

Total running time of the script: ( 1 minutes 38.118 seconds)

Gallery generated by Sphinx-Gallery

previous

Selecting random state for CSC

next

Extracting cross-frequency coupling waveforms from rodent LFP data