Stage 22: time-frequency power and inter-trial coherence

22-kit_time_frequency.py decomposes the stage-9 epochs of every recording and condition into time and frequency with Morlet wavelets (mne.Epochs.compute_tfr(method="morlet")), averages the power over trials and computes the inter-trial coherence (ITC). It answers what the evoked response cannot: oscillatory activity that rises or falls after an event without a fixed phase (alpha suppression over visual cortex when a stimulus appears, the beta rebound after a movement, gamma bursts) averages out in the time domain and is only visible in the trial-wise power.

python 22-kit_time_frequency.py -c pipeline_config_files/epilepsy_project_config.yml [--overwrite] [--no-plots]

Inputs: the stage-9 epochs files (*_ev<code>_epo.fif) and the condition names of events.conditions. Outputs, under derivatives/22-kit_time_frequency/sub-<id>/:

  • <recording>_cond-<name>_desc-power_tfr.npz: the power (sensors x frequencies x times), not baseline-corrected, with freqs, times, ch_names and nave; _desc-power_tfr.h5 next to it is MNE’s own AverageTFR file (written when h5io is installed); <recording>_desc-tfr-info.fif holds the measurement info;

  • <recording>_cond-<name>_desc-itc_tfr.npz (+ .h5): the inter-trial coherence;

  • <recording>_cond-<name>_desc-bandpower.tsv: time x bands, the baseline-corrected power averaged over sensors and over each band;

  • <recording>_desc-timefrequency.json: provenance (input identity, settings), the peaks per band and the ITC peak;

  • figures/*_desc-tfr.png (power, mean over sensors), *_desc-itc.png, *_desc-tfrbands.png (band time courses), *_desc-tfrtopo.png (one topomap per band at the band’s peak);

  • time_frequency_summary.csv: one row per recording and condition.

What it does

  1. Frequencies and wavelets. time_frequency.freqs gives n_freqs centre frequencies between fmin and fmax (log spacing by default). Each wavelet has cycles_per_hz x frequency cycles, at least min_cycles (n_cycles sets a fixed number instead). More cycles sharpen the frequency resolution and blur the time resolution. The Gaussian envelope of a wavelet at frequency f with n cycles has the width sigma_t = n / (2 pi f) s; a frequency is kept when 6 sigma_t (99.7 % of the envelope) fits into the epoch, the others are dropped and listed in freqs_dropped. With cycles_per_hz the envelope width is the same at every frequency (0.08 s for 0.5 cycles per Hz, so epochs of 0.5 s and longer keep every frequency; min_cycles widens the lowest ones). The first and last 3 sigma_t of the epoch mix with what lies beyond the epoch: this stretch is recorded as edge_s, hatched in the figures and excluded from the peak search. A response that falls in it needs longer epochs (epoching.tasks). MNE’s wavelets are 10 sigma_t long and compute_tfr refuses a wavelet longer than the signal, so an epoch shorter than that (0.8 s with 0.5 cycles per Hz) is mirrored at both ends up to the wavelet length before the convolution and cropped back afterwards (pad_s in the summary says by how much). Conditions with fewer than min_epochs trials (5 by default) are skipped (skipped_few_epochs).

  2. Power and ITC. epochs.compute_tfr(..., average=True, return_itc=True, decim=...): the power of every trial averaged (this keeps induced activity) and the ITC, the length of the mean of the unit phase vectors over trials (1 = the same phase in every trial, 0 = random; high where the response is evoked). The KIT gradiometers (typed mag) are used, bad channels excluded.

  3. Baseline correction for the figures and the band time courses (baseline.mode: logratio = log10 of power over the mean power in the baseline, 0 = no change, +0.3 = twice the baseline; percent, ratio, zscore, mean or null). baseline.window: null takes the epoch’s own baseline window from the stage-9 record. The saved power is not corrected, so any baseline can be applied later (power.apply_baseline((None, 0), mode="logratio")).

  4. Band time courses and peaks. For each band of time_frequency.bands the corrected power is averaged over sensors and over the band’s frequencies; the summary reports the band with the largest absolute change after the event (peak_band, peak_time_s, peak_value) and the time and frequency of the largest ITC (itc_peak, itc_peak_time_s, itc_peak_freq_hz).

Configuration

time_frequency:
  enabled: true
  tasks: []               # BIDS task labels; [] = every task with stage-9 epochs
  conditions: []          # condition names or codes; [] = all
  freqs: {fmin: 4, fmax: 40, n_freqs: 20, spacing: log}
  cycles_per_hz: 0.5      # n_cycles = cycles_per_hz x frequency ...
  min_cycles: 2           # ... but at least this many; or n_cycles: 7 for a fixed number
  n_cycles: null
  decim: 4                # keep every 4th sample of the result
  baseline: {mode: logratio, window: null}
  bands: {theta: [4, 8], alpha: [8, 13], beta: [15, 30]}
  itc: true
  max_epochs: 0           # 0 = all trials; a cap per condition for speed
  min_epochs: 5           # fewer trials: the condition is skipped
  save_plots: true
  n_jobs: 1

Reuse: a condition is recomputed when its epochs file changed (size and content digest), when the settings changed, or with --overwrite.

Reading the outputs in Python

import os
from pathlib import Path
from pipeline.mne_pipelines.kit_general_pipelines.results import Results
res = Results(bids_root=Path(os.environ["MEG_DATA"]) / "epilepsy-project")
rec = res.recording("0989", task="flippedsquares", run="01")
power = rec.tfr("flippedsquares_reversal_to_a")          # mne AverageTFR, not baseline-corrected
power.apply_baseline((None, 0), mode="logratio")
fig = power.plot(picks="mag", combine="mean", show=False)[0]
itc = rec.tfr("flippedsquares_reversal_to_a", kind="itc")
bands = rec.tfr_bands("flippedsquares_reversal_to_a")     # pandas: time + one column per band

Every figure of this stage in the report carries the same code under a “Reproduce this figure in Python” block (see Working with the results in Python (the results module and the report’s code blocks)).