| jupytext | kernelspec | ||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
Synopsis
Data: 100 s of local field potential data sampled at 1000 Hz.
Goal: Characterize the coupling between rhythms of different frequency.
Tools: Hilbert transform, analytic signal, instantaneous phase, cross-frequency coupling.
+++
+++
We begin this module with an "on-ramp" to analysis. The purpose of this on-ramp is to introduce you immediately to a core concept in this module: how to compute cross-frequency coupling (CFC) in Python. You may not understand all aspects of the program here, but that's not the point. Instead, the purpose of this on-ramp is to illustrate what can be done. Our advice is to simply run the code below and see what happens ...
from scipy.io import loadmat
from scipy import signal
from pylab import *
%matplotlib inline
rcParams['figure.figsize']=(12,3) # Change the default figure size
data = loadmat('matfiles/LFP-1.mat') # Load the LFP data,
t = data['t'][0] # ... extract t, the time variable,
LFP = data['LFP'][0] # ... and LFP, the voltage variable,
dt = t[1] - t[0] # Define the sampling interval,
fNQ = 1 / dt / 2 # ... and Nyquist frequency.
Wn = [5,7]; # Set the passband [5-7] Hz,
n = 100; # ... and filter order,
b = signal.firwin(n, Wn, nyq=fNQ, pass_zero=False, window='hamming');
Vlo = signal.filtfilt(b, 1, LFP); # ... and apply it to the data.
Wn = [80, 120]; # Set the passband [80-120] Hz,
n = 100; # ... and filter order,
b = signal.firwin(n, Wn, nyq=fNQ, pass_zero=False, window='hamming');
Vhi = signal.filtfilt(b, 1, LFP); # ... and apply it to the data.
phi = angle(signal.hilbert(Vlo)) # Compute phase of low-freq signal
amp = abs(signal.hilbert(Vhi)) # Compute amplitude of high-freq signal
p_bins = arange(-pi,pi,0.1) # To compute CFC, define phase bins,
a_mean = zeros(size(p_bins)-1) # ... variable to hold the amplitude,
p_mean = zeros(size(p_bins)-1) # ... and variable to hold the phase.
for k in range(size(p_bins)-1): # For each phase bin,
pL = p_bins[k] #... get lower phase limit,
pR = p_bins[k+1] #... get upper phase limit.
indices=(phi>=pL) & (phi<pR) #Find phases falling in this bin,
a_mean[k] = mean(amp[indices]) #... compute mean amplitude,
p_mean[k] = mean([pL, pR]) #... save center phase.
plot(p_mean, a_mean) #Plot the phase versus amplitude,
ylabel('High-frequency amplitude') #... with axes labeled.
xlabel('Low-frequency phase')
title('CFC');
Q: Try to read the code above. Can you see how it loads data, computes the phase and amplitude of the signals, and assess the cross-frequency coupling?
A: If you've never computed cross-frequency coupling before, that's an especially difficult question. Please continue on to learn this and more!
+++
from IPython.lib.display import YouTubeVideo
YouTubeVideo('Q-VQaY6iDMM')
In this module, we focus on local field potential (LFP) recordings. The LFP is a measure of local population neural activity, produced from small aggregates of neurons. In these data, we examine the association between rhythms of different frequencies.
In general, lower-frequency rhythms have been observed to engage larger brain areas and modulate spatially localized fast oscillations (see, for example). Of the different types of cross-frequency coupling (CFC) observced between brain rhythms, perhaps the most is phase-amplitude coupling (PAC), in which the phase of a low frequency rhythm modulates the amplitude (or power) of a high frequency rhythm. Cross-frequency coupling been observed in many brain regions, has been shown to change in time with task demands, and has been proposed to serve a functional role in working memory, neuronal computation, communication, and learning. Although the cellular and dynamic mechanisms of specific rhythms associated with CFC are relatively well understood, the mechanisms governing interactions between different frequency rhythms - and the appropriate techniques for measuring CFC - remain active research areas. Although we consider only a single electrode recording here, we note that these techniques can be extended to association measures between electrodes as well.
We are approached by a collaborator recording the local field potential (LFP) from rat hippocampus. She has implanted a small bundle of electrodes, which remain (chronically) implanted as the rat explores a large circular arena. She is interested in assessing the association between different frequency rhythms of the LFP, and more specifically whether an association between different frequency rhythms exists as the rat explores the arena. To address this question, she has provided us with 100 s of LFP data recorded during the experiment (i.e., while the rat spontaneously explored the arena).
+++
Our goal is to assess the associations between different frequency rhythms recorded in the LFP. To do so, we analyze the LFP data by computing the phase-amplitude coupling (CFC) of the time series. We construct two CFC measures that characterize how the phase of a low-frequency signal modulates the amplitude envelope of a high-frequency signal.
In this chapter, we develop two CFC measures, with a particular focus on phase-amplitude coupling (PAC). We introduce the concepts of the Hilbert transform, analytic signal, instantaneous phase, and amplitude envelope.
YouTubeVideo('mLglqyb55_Y')
Let's begin with visual inspection of the LFP data.
# Load the modules and set plot defaults
from scipy.io import loadmat
from pylab import *
%matplotlib inline
rcParams['figure.figsize']=(12,3) # Change the default figure size
data = loadmat('matfiles/LFP-1.mat') # Load the LFP data,
t = data['t'][0] # ... extract t, the time variable,
LFP = data['LFP'][0] # ... and LFP, the voltage variable,
plot(t, LFP) # ... and plot the trace,
xlabel('Time [s]') # ... with axes labeled.
ylabel('Voltage [mV]');
Next, let's take a closer look at an example 1 s interval of the data:
plot(t, LFP) # Plot the LFP data,
xlim([4, 5]) # ... restrict the x-axis to a 1 s interval,
ylim([-2, 2]) # ... and tighten the y-axis.
xlabel('Time [s]') # Label the axes
ylabel('Voltage [mV]')
show()
Visual inspection immediately suggests a dominant low-frequency rhythm interspersed with smaller-amplitude blasts of high-frequency activity.
Q. Approximate features of the rhythmic activity by visual inspection of the plot above. What is the frequency of the large-amplitude rhythm? Do you observe high-frequency activity? If so, where in time, and at what approximate frequency? What is the sampling frequency of these data? If you were to compute the spectrum of the entire dataset (100 s of LFP), what would be the Nyquist frequency and the frequency resolution? Hint: Consider the times near 4.35 s and 4.5 s. Do you see the transient fast oscillations?
+++
Visual inspection of the LFP data suggests that multiple rhythms appear. To further characterize this observation, we compute the spectrum of the LFP data. We analyze the entire 100 s of data and compute the spectrum with a Hanning taper.
dt = t[1] - t[0] # Define the sampling interval,
T = t[-1] # ... the duration of the data,
N = len(LFP) # ... and the no. of data points
x = hanning(N) * LFP # Multiply data by a Hanning taper
xf = rfft(x - x.mean()) # Compute Fourier transform
Sxx = 2*dt**2/T * (xf*conj(xf)) # Compute the spectrum
Sxx = real(Sxx) # Ignore complex components
df = 1 / T # Define frequency resolution,
fNQ = 1 / dt / 2 # ... and Nyquist frequency.
faxis = arange(0, fNQ + df, df) # Construct freq. axis
plot(faxis, 10 * log10(Sxx)) # Plot spectrum vs freq.
xlim([0, 200]) # Set freq. range,
ylim([-80, 0]) # ... and decibel range
xlabel('Frequency [Hz]') # Label the axes
ylabel('Power [mV$^2$/Hz]');
The resulting spectrum reveals two intervals of increased power spectral density. The lowest-frequency peak at 6 Hz is also the largest and corresponds to the slow rhythm we observe dominating the signal through visual inspection. Remember a plot of that signal:
At higher frequencies, we find an additional broadband peak at approximately 80–120 Hz. These spectral results support our initial visual inspection of the signal; there exist both low- and high-frequency activities in the LFP data. We now consider the primary question of interest: Do these different frequency rhythms exhibit associations?
+++
YouTubeVideo('JjOcJy4dVzE')
To assess whether different frequency rhythms interact in the LFP recording, we implement a measure to calculate the CFC. The idea of CFC analysis, in our case, is to determine whether a relation exists between the phase of a low-frequency signal and the envelope or amplitude of a high-frequency signal. In general, computing CFC involves three steps. Each step contains important questions and encompasses entire fields of study. Our goal in this section is to move quickly forward and produce a procedure we can employ, investigate, and criticize. Continued study of CFC - and the associated nuances of each step - is an active area of ongoing research.
-
Filter the data into high- and low-frequency bands.
-
Extract the amplitude and phase from the filtered signals.
-
Determine if the phase and amplitude are related.
+++
YouTubeVideo('WL_nFRBHqLU')
The first step in the CFC analysis is to filter the data into two frequency bands of interest. The choice is not arbitrary: the separate frequency bands are motivated by initial spectral analysis of the LFP data. In this case, we choose the low-frequency band as 5–7 Hz, consistent with the largest peak in the spectrum, and the high-frequency band as 80–120 Hz, consistent with the second-largest broadband peak. To consider alternative frequency bands, the same analysis steps would apply.
+++
There are many options to perform the filtering. To do so requires us to design a filter that ideally extracts the frequency bands of interest without distorting the results. Here, we apply a finite impulse response (FIR) filter. We start by extracting the low-frequency band:
from scipy import signal
Wn = [5,7]; # Set the passband [5-7] Hz,
n = 100; # ... and filter order,
# ... build the bandpass filter,
b = signal.firwin(n, Wn, nyq=fNQ, pass_zero=False, window='hamming');
Vlo = signal.filtfilt(b, 1, LFP); # ... and apply it to the data.
Next, we extract the high-frequency band:
Wn = [80, 120]; # Set the passband [80-120] Hz,
n = 100; # ... and filter order,
# ... build the bandpass filter,
b = signal.firwin(n, Wn, nyq=fNQ, pass_zero=False, window='hamming');
Vhi = signal.filtfilt(b, 1, LFP); # ... and apply it to the data.
For each frequency band, we specify a frequency interval of interest by defining the low- and high-cutoff frequencies in the variable Wn. In this way, we specify the passband of the filter. We then set the filter order (n) and design the filter using the Python function signal.firwin from the SciPy package. Finally, we apply the filter using the Python function signal.filtfilt, which performs zero-phase filtering by applying the filter in both the forward and reverse directions. We note that, in this case, the filtering procedure is nearly the same in both frequency bands; the only change is the specification of the frequency interval of interest.
To understand the impact of this filtering operation on the LFP, let’s plot the results. More specifically, let's plot the original signal, and the signal filtered in the low- and high-frequency bands, for an example 2 s interval of time:
figure(figsize=(14, 4)) # Create a figure with a specific size.
plot(t, LFP) # Plot the original data vs time.
plot(t, Vlo) # Plot the low-frequency filtered data vs time.
plot(t, Vhi) # Plot the high-frequency filtered data vs time.
xlabel('Time [s]')
xlim([24, 26]); # Choose a 2 s interval to examine
ylim([-2, 2]);
legend(['LFP', 'Vlo', 'Vhi']); # Add a legend.
As expected, the low-frequency band (orange) captures the large-amplitude rhythm dominating the LFP signal, while the higher-frequency band (green) isolates the brief bursts of faster activity.
+++
YouTubeVideo('QiRuBvbCQt4')
The next step in the CFC procedure is to extract the phase of the low-frequency signal and the amplitude envelope (or simply, amplitude) of the high-frequency signal. To compute a particular type of CFC, namely phase-amplitude coupling, we then compare these two signals, i.e., we compare the phase of the low-frequency activity and the amplitude envelope of the high-frequency activity. How do we actually extract the phase and amplitude signals from the data? There are a variety of options to do so, and we choose here to employ the analytic signal approach, which allows us to estimate the instantaneous phase and amplitude envelope of the LFP.
YouTubeVideo('Ir8Gf550S3o')
The first step in computing the analytic signal is to compute the Hilbert transform. We begin with some notation. Define
It’s perhaps more intuitive to consider the effect of the Hilbert Transform on the frequencies
To summarize: The Hilbert transform
The Hilbert Transform can also be described in the time domain, although its representation is hardly intuitive (see this Appendix for more details). Then the analytic signal
The effect of the Hilbert transform is to remove negative frequencies from
+++
YouTubeVideo('-CjnFEOopfw')
Let
where to simplify notation we have defined
The real variable
For real signals, which include nearly all recordings of brain activity, the negative frequency component is redundant, and we usually ignore it. However, the negative frequency component still remains. By definition, the effect of the Hilbert transform is to induce a phase shift. For positive frequencies, the phase shift is
+++
Q: Why does a phase shift of
A: Consider the complex exponential
Notice the result simplifies to the original complex exponential
+++
Q. Can you show that a
+++
This analysis shows that we can represent the Hilbert transform of
Therefore, the Hilbert transform of
In this equation, we multiply the positive frequency part of
Q. Can you perform this simplification? In other words, can you show that
YouTubeVideo('e4kECy8KP-4')
The Hilbert transform of
YouTubeVideo('emsU97RcFjI')
We are now ready to define the analytic signal
Notice that this analytic signal
compared to only one complex exponential required to express the corresponding analytic signal
Because
Q. The phase of a complex quantity is the angle with respect to the real axis in the complex plain. What is the angle of
+++
Having developed some understanding of the Hilbert Transform, let’s now return to the LFP data of interest here. It’s relatively easy to compute the analytic signal and extract the phase and amplitude in Python:
phi = angle(signal.hilbert(Vlo)) # Compute phase of low-freq signal
amp = abs(signal.hilbert(Vhi)) # Compute amplitude of high-freq signal
These operations require just two lines of code. But beware the following.
Alert! The command hilbert(x) returns the analytic signal of
To extract the phase, we apply the function angle to the analytic signal of the data filtered in the low-frequency band (variable Vlo). We then extract the amplitude of the analytic signal of the data filtered in the high-frequency band (variable Vhi) by computing the absolute value (function abs).
To summarize, in this step we apply the Hilbert transform to create the analytic signal and get the phase or amplitude of the bandpass-filtered data.
+++
YouTubeVideo('fiL9UNbLPA8')
As with the previous steps, we have at our disposal a variety of procedures to assess the relation between the phase (of the low-frequency signal) and amplitude (of the high-frequency signal). We do so here in one way, by computing the phase-amplitude plot.
+++
To start, define the two-column phase-amplitude vector,
where
We now use this two-column vector to assess whether the phase and amplitude envelope are related. Let's begin by plotting the average amplitude versus phase. We divide the phase interval into bins of size 0.1 beginning at
p_bins = arange(-pi, pi, 0.1)
a_mean = zeros(size(p_bins)-1)
p_mean = zeros(size(p_bins)-1)
for k in range(size(p_bins)-1): #For each phase bin,
pL = p_bins[k] #... lower phase limit,
pR = p_bins[k+1] #... upper phase limit.
indices=(phi>=pL) & (phi<pR) #Find phases falling in bin,
a_mean[k] = mean(amp[indices]) #... compute mean amplitude,
p_mean[k] = mean([pL, pR]) #... save center phase.
plot(p_mean, a_mean) #Plot the phase versus amplitude,
ylabel('High-frequency amplitude') #... with axes labeled.
xlabel('Low-frequency phase');
Q. Consider this plot of the average amplitude versus phase. Does this result suggest CFC occurs in these data?
A. The plot shows the phase bins versus the mean amplitude in each bin. Visual inspection of this phase-amplitude plot suggests that the amplitude of the high-frequency signal depends on the phase of the low-frequency signal. In particular, we note that when the phase is near a value of 2 radians, the amplitude tends to be large, while at other phases the amplitude is smaller. This conclusion suggests that CFC does occur in the data; the high-frequency amplitude depends on the low-frequency phase.
+++
Q. If no CFC occurred in the data, what would you expect to find in the plot of average amplitude versus phase?
+++
As a basic statistic to capture the extent of this relation, we compute the difference between the maximum and minimum of the average amplitude envelope over phases. Let's assign this difference the label
h = max(a_mean)-min(a_mean)
print(h)
We find a value of
For each surrogate phase-amplitude vector, we compute the statistic
n_surrogates = 1000; #Define no. of surrogates.
hS = zeros(n_surrogates) #Vector to hold h results.
for ns in range(n_surrogates): #For each surrogate,
ampS = amp[randint(0,N,N)] #Resample amplitude,
p_bins = arange(-pi, pi, 0.1) #Define the phase bins
a_mean = zeros(size(p_bins)-1) #Vector for average amps.
p_mean = zeros(size(p_bins)-1) #Vector for phase bins.
for k in range(size(p_bins)-1):
pL = p_bins[k] #... lower phase limit,
pR = p_bins[k+1] #... upper phase limit.
indices=(phi>=pL) & (phi<pR) #Find phases falling in bin,
a_mean[k] = mean(ampS[indices]) #... compute mean amplitude,
p_mean[k] = mean([pL, pR]) #... save center phase.
hS[ns] = max(a_mean)-min(a_mean) # Store surrogate h.
In this code, we first define the number of surrogates (variable n_surrogates) and a vector to store the statistic hS). Then, for each surrogate, we use the Python function randint to randomly permute the amplitude vector without replacement. We then use this permuted amplitude vector (variable ampS) and the original phase vector (variable phi) to compute h for the original (unpermuted) data.
Let's now view the results of this resampling procedure by creating a histogram of the variable hS, and compare this distribution to the value of h we computed earlier.
counts, _, _ = hist(hS, label='surrogates') # Plot the histogram of hS, and save the bin counts.
vlines(h, 0, max(counts), colors='red', label='h', lw=2) # Plot the observed h,
legend(); # ... include a legend.
The value of h) lies far outside the distribution of surrogate hS). To compute a
p = sum([s > h for s in hS]) / len(hS)
print(p)
We find a
from IPython.lib.display import YouTubeVideo
YouTubeVideo('NQUfELSZ0Cc')
from IPython.lib.display import YouTubeVideo
YouTubeVideo('NQUfELSZ0Cc')
In this module, we considered techniques to characterize cross-frequency coupling (CFC), associations between rhythmic activity observed in two different frequency bands. To do so, we introduced the Hilbert transform, which can be used to compute the instantaneous phase and amplitude of a signal. We focused on characterizing the relation between the phase of low-frequency band activity (5-7 Hz) and the amplitude of high-frequency band activity (100-140 Hz). To do this, we computed the average amplitude at each phase and determined the extent of the variability (or wiggliness).
For the LFP data of interest here, we found evidence of CFC between the two frequency bands. Importantly, these results also appear consistent with visual inspection of the unfiltered data. Careful inspection of the example voltage traces suggests that CFC does in fact occur in these data. In general, such strong CFC, visible to the naked eye in the unprocessed LFP data, is unusual. Instead, data analysis techniques are required to detect features not obvious through visual inspection alone. Developing techniques to assess CFC and understanding the biophysical mechanisms of CFC and implications for function, remain active research areas.
In developing these approaches, we utilized expertise and procedures developed in other modules. In particular, we relied on the notions of frequency and power, amplitude and phase, filtering, and resampling. Such a multifaceted approach is typical in analysis of neural data, where we leverage the skills developed in analyzing other datasets.
+++
If you enjoy Case-Studies-Python, and would like to share your enjoyment with us, sponsor our coffee consuption here.
+++
We have presented the Hilbert Transform in the frequency domain: it produces a quarter cycle phase shift. It's reasonable to consider as well the time domain representation of the Hilbert Transform. To do so, let's write the Hilbert Transform as
where we have written
In the frequency domain, we perform the Hilbert Transform by multiplying the signal
To convert the Hilbert transform in the frequency domain to the Hilbert Transform in the time domain, we take the inverse Fourier transform. Looking up the inverse Fourier transform of
Let's represent the inverse Fourier transform of
Now, let's make use of an important fact. Multiplication of two quantities in the frequency domain corresponds to convolution of these two quantities in the time domain (see notebook 04). The convolution of two signals
So, in the time domain, the Hilbert transform becomes
This time domain representation of the HIlbert Transform is equivalent to the frequency domain representation. However, the time domain representation is much less intuitive. Compare the previous equation to the statement, "The Hilbert Transform is a 90-degree phase shift in the frequency domain." The latter, we propose, is much more intuitive.