Oscillation Presence

Topic: neural oscillations are variably present, and therefore have to be explicitly detected.

Issue

Neural activity contains both periodic (oscillatory) and aperiodic activity, which contributes power across all frequencies.

Neural oscillations (periodic activity) are variable in their occurrence. Due to aperiodic activity, there will always be power at any given frequency.

Solution

Investigations of neural oscillations should include a detection step, in which the presence of oscillatory activity is explicitly validated.

This can be done be with spectral parameterization, as the presence of a frequency specific oscillation can be detected as a peak in the power spectrum.

Settings

import seaborn as sns
sns.set_context('poster')
# Set random seed
set_random_seed(808)
# Define general simulation settings
n_seconds = 20
fs = 1000

times = create_times(n_seconds, fs)
# Define parameters for the simulations
exp = -1.5
cf = 10
ap_filt = (1.5, 75)
# Set frequency ranges of interest
plt_range = [2, 30]
filt_xlim = [5, 8]
# Plot settings
plt_kwargs = {'xlabel' : '', 'ylabel' : '', 'lw' : 3}
shade_colors = [BAND_COLORS[band] for band in BANDS.labels]
# Set whether to save out figures
SAVE_FIG = False

Aperiodic Signal Have Power Across all Frequencies

For an initial example, let’s consider the dirac delta - a single value of 1 in a series of zeros.

# Simulate a delta function
n_points = len(times)
dirac_sig = np.zeros([n_points])
dirac_sig[round(n_points/2)] = 1
# Plot the time series of the delta signal
plot_time_series(times, dirac_sig, lw=4)
plt.axis('off');
savefig(SAVE_FIG, '01-dirac_ts')
../_images/01-OscillationPresence_13_0.png
# Compute a power spectrum of the Dirac delta
freqs, powers = compute_spectrum(dirac_sig, fs)
plot_spectrum(*trim_spectrum(freqs, powers, [1.5, 30]), 
              color='black', lw=6, log_powers=False)

plt.xlim(plt_range)
style_psd(plt.gca())
savefig(SAVE_FIG, '01-dirac_psd')
../_images/01-OscillationPresence_14_0.png

Rhythmic signals have power at specific frequencies

Next, lets examine a pure sinusoid.

# Simulate a sinusoid
osc_sig = sim_oscillation(n_seconds, fs, 10)
# Plot the time series of the sinusoidal signal
plot_time_series(times, osc_sig, lw=3, xlim=[2, 3])
plt.axis('off');
savefig(SAVE_FIG, '01-sine-ts')
../_images/01-OscillationPresence_17_0.png
# Compute a power spectrum of the Dirac delta
freqs, powers = compute_spectrum(osc_sig, fs, nperseg=2*fs)
plot_spectrum(*trim_spectrum(freqs, powers, [1, 50]), 
              color='black', lw=6, log_powers=False)

plt.xlim(plt_range)
style_psd(plt.gca())

savefig(SAVE_FIG, '01-sine-psd')
../_images/01-OscillationPresence_18_0.png

Simulate Aperiodic Time Series

First, let’s simulate an example time series.

For this example, we want to emphasize that measures may not necessarily reflect oscillatory activity.

To do so, we will simulate a signal with neurally plausible aperiodic activity, but with no oscillations.

In this scenario, we can examine measurements that reflect only the aperiodic activity.

# Simulate aperiodic time series
ap_sig = sim_powerlaw(n_seconds, fs, exponent=exp, f_range=ap_filt)
# Plot the simulated time series
_, ax = plt.subplots(figsize=(7, 2))
plot_time_series(times, ap_sig, xlim=(10, 13), alpha=0.75, ax=ax)
plt.axis('off')
savefig(SAVE_FIG, '01-ts')
../_images/01-OscillationPresence_21_0.png
# Zoom in and plot a segment of the simulated signal
st, en = 10.2, 11
_, ax = plt.subplots(figsize=(8, 2.5))
plot_time_series(times, ap_sig, ax=ax, xlim=[st, en], **plt_kwargs)
ax.set_xticks([]); ax.set_yticks([]);
savefig(SAVE_FIG, '01-ts-segment')
../_images/01-OscillationPresence_22_0.png

Compute Power Spectrum

Next, we can compute the power spectrum of our measured signal.

Since we simulated a signal with aperiodic activity, note that we expect the power spectrum to reflect this 1/f property of the data.

# Compute power spectrum of the simulated signal
freqs, powers = trim_spectrum(*compute_spectrum(ap_sig, fs, nperseg=fs), plt_range)
# Plot the power spectrum of the simulated signal
plot_spectrum_shading(freqs, powers, lw=5, color='black', log_powers=True,
                      shades=BANDS.definitions, shade_colors=shade_colors)
plt.xlim(plt_range)
style_psd(plt.gca())
savefig(SAVE_FIG, '01-ap_psd_col')
../_images/01-OscillationPresence_25_0.png

In the spectrum above we do see the expected 1/f nature of the power spectrum.

Notably, due to the 1/f there is power across all frequencies, despite there being no oscillatory power.

Examine Filtered Ranges

Next, let’s examine what the data looks like when filtered in band-specific frequency ranges.

# Create a plot of the data filtered into different frequency ranges
_, axes = plt.subplots(BANDS.n_bands, 1, figsize=(8, 9))
for ax, color, (label, f_range) in zip(axes, shade_colors, BANDS):

    band_sig = filter_signal(ap_sig, fs, 'bandpass', f_range)
    plot_time_series(times, band_sig, ax=ax,
                     xlim=filt_xlim, ylim=(-1.25, 1.25), **plt_kwargs)
    ax.axvspan(filt_xlim[0], filt_xlim[1], alpha=0.2, color=color)
    ax.set_xticks([]); ax.set_yticks([]); 
    
plt.subplots_adjust(hspace=0.05)
savefig(SAVE_FIG, '01-ap_ts_bands')
../_images/01-OscillationPresence_28_0.png

Note in the above filtered data, that there is activity in all bands, and it even seems to have interesting dynamics.

This is because the aperiodic activity has activity at all frequencies, that can be extracted with filters.

In addition, the seeming sinusoidality of the filtered traces stems from the filters themselves.

However, these traces do not reflect true oscillatory activity, and interpreting it as such would be wrong.

Compare to a signal with an oscillation

# Simulate a time series with an oscillation
osc_sig = sim_oscillation(n_seconds, fs, freq=cf, variance=0.5)

# Combine the oscillation with the aperiodic signal from before
sig_comb = ap_sig + osc_sig
# Plot the simulated time series
plot_time_series(times, sig_comb, xlim=(10, 13), lw=2.5, alpha=0.75)
plt.axis('off')
savefig(SAVE_FIG, '01-comb_ts')
../_images/01-OscillationPresence_32_0.png
# Compute power spectrum of the simulated signal
freqs, powers = trim_spectrum(*compute_spectrum(sig_comb, fs, nperseg=1.5*fs), plt_range)
# Plot the power spectrum of the simulated signal
plot_spectrum_shading(freqs, powers, lw=5, color='black', log_powers=True,
                      shades=BANDS.definitions, shade_colors=shade_colors)
plt.xlim(plt_range)
style_psd(plt.gca())
savefig(SAVE_FIG, '01-comb_psd_col')
../_images/01-OscillationPresence_34_0.png

Peak detection

In order to detect an oscillation, we can examine the the power spectrum for a peak of frequency-specific activity.

Here, we will use spectral parameterization to detect peaks in the power spectrum.

# Initialize an object for spectral parameterization
fm = FOOOF(verbose=False)

# Fit the spectral parameterization model
fm.fit(freqs, powers)
# Set the plot color for the shaded peak to be the alpha color
from fooof.plts.settings import PLT_COLORS
PLT_COLORS['periodic'] = shade_colors[2]
# Plot the power spectrum, with the model, highlighting the detected peak
plot_annotated_model(fm, annotate_peaks=False, annotate_aperiodic=False)
plt.gca().get_legend().remove()
plt.xlim(plt_range)
style_psd(plt.gca())
savefig(SAVE_FIG, '01-comb_psd_peak')
../_images/01-OscillationPresence_38_0.png

Having detected a peak, we can restrict our analyses to frequency bands with verified oscillations.

# Plot the filtered range for the band with a detected peak
_, ax = plt.subplots(figsize=(8, 3))
band_sig = filter_signal(sig_comb, fs, 'bandpass', BANDS.alpha)
plot_time_series(times, band_sig, ax=ax,
                 xlim=filt_xlim, ylim=(-1.5, 1.5), **plt_kwargs)
ax.set_xticks([]); ax.set_yticks([]); 
ax.axvspan(filt_xlim[0], filt_xlim[1], alpha=0.2, color=shade_colors[2])
savefig(SAVE_FIG, '01-alpha_ts_bands')
../_images/01-OscillationPresence_40_0.png

Conclusion

Neural often, but not always, contains periodic activity.

Since there is a lot of other activity in neural recordings, measures of particular frequency ranges will always return some activity, even if it is not oscillatory.

An explicit detection step should be done, in which the presence of oscillatory activity is validated.