.ipynb .pdf

Binder

Contents

Waveform Shape

Neural oscillations vary in their waveform shape, which can impact their measurements and interpretations.

Issue

Neural oscillations are often non-sinusoidal, and exhibit variability in their waveform shape.

This can cause issues with analysis methods that assume sinusoidal bases.

This includes confounded power measurements (due to harmonics), spurious phase-amplitude coupling (PAC), and biased filter outputs.

Solution

Waveform shape can and should be explicitly checked and measured.

Related Work

Core citations / links:

• When brain rhythms aren’t ‘rhythmic’: implication for their mechanisms and meaning, Jones, 2016

• Brain Oscillations and the Importance of Waveform Shape, Cole & Voytek, 2017

Related Discussion / Examples

This example uses the ByCycle module for measuring cycle-by-cycle variation.

• Tutorials on the ByCycle approach are available here

• A relevant example on the interaction of shape and PAC is here

import numpy as np
import matplotlib.pyplot as plt

from fooof.plts.spectra import plot_spectra, plot_spectrum_shading
from fooof.utils import trim_spectrum

from neurodsp.filt import filter_signal
from neurodsp.spectral import compute_spectrum
from neurodsp.plts.time_series import plot_time_series
from neurodsp.utils import create_times, set_random_seed
from neurodsp.sim import sim_oscillation, sim_combined, sim_powerlaw

# Import bycycle module, which is used for this example
from bycycle.features import compute_features

# Add path and import custom code for the project
import sys
sys.path.append('code')
from utils import mu_wave, compute_pac
from plts import (plot_waveforms, plot_power_by_shape, plot_harmonic_power, 
                  plot_pac, plot_legend, savefig)
from settings import ALPHA_COLOR, ALPHA_RANGE

Settings

import seaborn as sns
sns.set_context('talk')

# Set random seed
set_random_seed(808)

# Set the average function to use
avg_func = np.nanmean

# Define general simulation settings
n_seconds = 25
fs = 1000

times = create_times(n_seconds, fs)

# Define parameters for the simulations
exp = -1.5
ap_filt = (2, 150)
cf = 10
rdsym = 0.05

# Define ranges of interest
psd_range = (2, 50)
beta_range = (15, 35)

# Define range of values to use for oscillation asymmetry
rdsyms = [0.50, 0.625, 0.75, 0.875]

# Collect parameters and set up simulations
comps = {'sim_powerlaw' : {'exponent' : exp, 'f_range' : ap_filt},
         'sim_oscillation' : {'freq' : cf, 'cycle' : 'asine', 'rdsym' : rdsym}}

# Define relative power of the signal components
comp_vars = [1, 0.8]

# Plot settings
plt_kwargs = {'xlabel' : '', 'ylabel' : '', 'alpha' : [0.75, 0.75]}

# Define plot colors
cmap = [plt.cm.viridis(i) for i in np.linspace(0, 1, len(rdsyms) + 1)]
cmap = cmap[::-1]

mu_color = "k"
sinus_color = "r"
harmonic_color = "b"
colors = [mu_color, sinus_color, harmonic_color]

# Burst detection settings for ByCycle
burst_detection_kwargs = {
    "amplitude_fraction_threshold": 0.0,
    "amplitude_consistency_threshold": 0.0,
    "period_consistency_threshold": 0.0,
    "monotonicity_threshold": 0.0}

# Set whether to save out figures
SAVE_FIG = False

Simulate Time Series

First, let’s visualize an asymmetric oscillation.

# Simulate an asymmetric oscillation
osc = sim_oscillation(n_seconds, fs, cf, cycle='asine', rdsym=rdsym)

# Visualize the asymmetric oscillation
plot_time_series(times, osc, xlim=[5, 6], **plt_kwargs)
plt.axis('off');
savefig(SAVE_FIG, '05-shape_osc')

# Simulate the combined time series
ap = sim_powerlaw(n_seconds, fs, exponent=exp, f_range=ap_filt, variance=0.1)

# Create the combiend signal
sig = osc + ap

# Plot the simulated time series
plot_time_series(times, sig, xlim=[5, 6.5], **plt_kwargs)
plt.axis('off')
savefig(SAVE_FIG, '05-ts_comb')

Computing Shape Features

To measure shape features, we will use the bycycle toolbox.

# Compute shape features of the simulated signal
df_shape = compute_features(sig, fs, f_range=(cf - 2, cf +2 ),
                            burst_detection_kwargs=burst_detection_kwargs)

# Check out some example features computed for the signal
df_shape.head()

	sample_peak	sample_zerox_decay	sample_zerox_rise	sample_last_trough	sample_next_trough	period	time_peak	time_trough	volt_peak	volt_trough	...	volt_rise	volt_amp	time_rdsym	time_ptsym	band_amp	amp_fraction	amp_consistency	period_consistency	monotonicity	is_burst
0	306	350	300	294	393	99	50	54	1.624848	-1.361975	...	2.986823	2.971948	0.121212	0.480769	1.259020	0.248980	NaN	NaN	0.866279	False
1	414	449	400	393	481	88	49	50	1.904694	-1.332224	...	3.236918	3.251017	0.238636	0.494949	1.338579	0.624490	0.888784	0.814815	0.801515	True
2	503	552	499	481	589	108	53	50	1.541561	-1.360423	...	2.901983	3.150959	0.203704	0.514563	1.184661	0.485714	0.853541	0.814815	0.686275	True
3	608	666	599	589	697	108	67	47	1.387864	-1.858374	...	3.246239	3.256740	0.175926	0.587719	1.121370	0.640816	0.869941	0.861111	0.741793	True
4	711	753	700	697	790	93	53	34	1.876328	-1.879377	...	3.755705	3.816923	0.150538	0.609195	1.181158	0.983673	0.869941	0.861111	0.814103	True

5 rows × 23 columns

Filter Time Series

Next, let’s explore filtering our simulated signal.

With the filtering we want to pay particular attention to how narrowband filters distort the shape of asymmetric oscillations.

# Filter the signal, including a broadband filter, and the alpha range
sig_filt_al = filter_signal(sig, fs, 'bandpass', ALPHA_RANGE)

# Plot the alpha filtered signal, compared to the original signal
plot_time_series(times, [sig, sig_filt_al], xlim=[5, 6.5], **plt_kwargs)
plt.axis('off')
savefig(SAVE_FIG, '05-ts_filt')

# Compute and plot the power spectrum of the signal
freqs1, powers1 = trim_spectrum(*compute_spectrum(sig, fs, nperseg=fs), (5, 35))
plot_spectrum_shading(freqs1, powers1, ALPHA_RANGE, color='k',
                      log_freqs=False, log_powers=True,
                      lw=4, shade_colors=ALPHA_COLOR)
savefig(SAVE_FIG, '05-psd')

Simulate Across Waveform Shapes

In the next set of simulations, we will create time series across variations of waveform shape.

# Simulate a collections of oscillations across variable waveforms shapes
oscs = []
for rdsym in rdsyms:
    oscs.append(sim_oscillation(n_seconds, fs, cf, cycle="asine", rdsym=rdsym))

# Plot the waveforms, comparing to filtered versions, to see how peak positions change
plot_waveforms(times, oscs, fs, cf, burst_detection_kwargs, cmap)
savefig(True, '05-ts_waveforms')

Compute Power Spectra

Next, let’s compute the power spectrum of the simulated signals.

In the power spectrum, we will be paying attention to harmonics that arise due to the asymmetric oscillation.

Note how, in the following plot, there are are clear harmonics in the beta range.

# Create an aperiodic component to add to signals
ap = sim_powerlaw(n_seconds, fs, exponent=-2.5)

# Compute power spectra for each signal, adding an aperiodic component
all_pows = []
for osc in oscs:
    cur_freqs, cur_pows = compute_spectrum(osc + 20 * ap, fs, nperseg=1.5 * fs)
    all_pows.append(cur_pows)

# Restrict power specra to range of interest
freqs, psd = trim_spectrum(cur_freqs, np.array(all_pows), psd_range)

# Plot spectra
_, ax_psd = plt.subplots(figsize=(5, 3))
for ind in range(4):
    ax_psd.semilogy(freqs, psd[ind].T, color=cmap[ind], lw=2)
ax_psd.set(xlabel="frequency [Hz]", xlim=(5, 35), yticks=[], ylabel="log PSD")
sns.despine(ax=ax_psd)
savefig(SAVE_FIG, '05-psd_waveforms')

Compare across asymmetry values

Next, let’s extend the analysis, and simulate signals across different oscillation asymmetries.

# Settings for current simulations
n_trials = 50
rdsyms_all = np.arange(0.5, 1, 0.025)
beta_pows = np.zeros((len(rdsyms_all), n_trials))

# Simulate for a number of trials
for i_rdsym, cur_rdsym in enumerate(rdsyms_all):
    for i_trial in range(n_trials):

        # Create the signal
        comps = {"sim_powerlaw": {"exponent": exp, "f_range": ap_filt},
                 "sim_oscillation": {"freq": cf, "cycle": "asine", "rdsym": cur_rdsym}}
        cur_sig = sim_combined(n_seconds, fs, comps, comp_vars)

        # Compute the spectrum and collect the measured beta power
        cur_freqs, cur_pows = compute_spectrum(cur_sig, fs, nperseg=fs)
        _, beta_pow = trim_spectrum(cur_freqs, cur_pows, beta_range)
        beta_pow = np.mean(beta_pow)
        beta_pows[i_rdsym, i_trial] = beta_pow

mean_beta = np.mean(beta_pows, axis=1)
std_beta = np.std(beta_pows, axis=1)

# Plot the measured beta power across waveform shape
plot_power_by_shape(rdsyms_all, mean_beta, std_beta, rdsyms, cmap)
savefig(SAVE_FIG, '05-power_by_shape')

In the above, we can see that the the measured beta power tracks the waveform shape of the simulated oscillation.

Keep in mind this is relating measured beta power to the shape of an alpha oscillation.

In this simulation, there is no actual beta power. This shows spurious power, due to the harmonics.

Compute Phase-Amplitude Coupling

Finally, let’s examine phase-amplitude coupling.

Phase-amplitude coupling is an analysis to investigate coupling between multiple, distinct oscillations.

Note that in our example signal, there is only one oscillation, and so there can be no true coupling between multiple periodic components.

Despite this, non-sinusoidal rhythms can induce artifactual phase-amplitude coupling.

# Define time vector for these simulations
time = np.arange(0, 2, 0.001)

# Simulate signals
signal_alpha, signal_beta = mu_wave(time, 0, cf, comb=False)
signal_mu = signal_alpha + signal_beta

# Filter signals to ranges of interest
filt_range = (15, 25)
signal_mu_filt = filter_signal(signal_mu, fs, "bandpass",
                               filt_range, remove_edges=False)
signal_beta_filt = filter_signal(signal_beta, fs, "bandpass",
                                 filt_range, remove_edges=False)

# Plot the time series and filtered traces
plot_harmonic_power(time, signal_mu, signal_beta, signal_mu_filt, colors)
savefig(SAVE_FIG, '05-pac_harmonic_power')

# Compute power spectra for different signal components
freqs, spec_mu = compute_spectrum(signal_mu, fs, f_range=psd_range)
freqs, spec_beta = compute_spectrum(signal_beta, fs, f_range=psd_range)
freqs, spec_mu_filt = compute_spectrum(signal_mu_filt, fs, f_range=psd_range)

# Plot a spectrum of the different signal components
_, axins = plt.subplots(figsize=(4.5, 3))
axins.semilogy(freqs, spec_mu, color=mu_color)
axins.semilogy(freqs, spec_beta, color=sinus_color)
axins.semilogy(freqs, spec_mu_filt, color=harmonic_color)
axins.set(xlim=(0, psd_range[1]), xticks=[], yticks=[]);
savefig(SAVE_FIG, '05-pac_psd_inset')

# Compute phase amplitude coupling
bins, pac = compute_pac(signal_mu_filt, signal_beta_filt, signal_alpha)

# Plot the calculate phase-amplitude coupling
plot_pac(bins, pac, [harmonic_color, sinus_color])
savefig(SAVE_FIG, '05-pac_coupling')

Conclusion

Many of the methods that are commonly used assume sinusoidality.

This can be an issue, as the observed oscillations in neural data are often non-sinusoidal.

As we’ve explored here, non-sinusoidal oscillations can induce spurious or biased measures of power, phase-amplitude coupling, and filter outputs.

These issues can be addressed by explicitly measuring waveform shape, both as a measure of interest itself, and to check whether waveform shape measures may explain any other measured changes.

previous

Temporal Fluctuations

next

Separation of Rhythms