Differences

This shows you the differences between two versions of the page.

--- analysis:course:week7 [2013/10/27 22:34]
mvdm [Application to real data]
+++ analysis:course:week7 [2018/07/07 10:19] (current)
@@ Line 1: / Line 1: @@
-===== Time-frequency analysis II: second-order relationships, phase-amplitude relationships =====
+===== Time-frequency analysis II: cross-frequency coupling =====
 Goals:
@@ Line 5: / Line 5: @@
   * Use the spectrogram as a basis to reveal correlations in power across frequencies
   * Understand how to compute and interpret autocorrelation and crosscorrelation functions
-  * Apply different methods to obtain signal phase for a given frequency
+  * Apply the Hilbert transform to obtain signal phase over time for a given frequency
-  * Relate phase of one frequency to power at another
+  * Relate phase of one frequency to power at another by using ''histc()''
 Deliverables:
-  * A function averageBin() that computes the average of signal B based on bins defined over signal A
+  * A function averageXbyYbin() that computes the average of signal B based on bins defined over signal A
-Resources:
+Assignment:
-  * 1
+  * Test for coupling between delta phase and low-gamma amplitude in a ventral striatal LFP
-  * 2
+Resources:
+  * (background reading, a short review on cross-frequency coupling) Jensen and Colgin, TICS 2007
+  * (optional, a nice application of the tools covered here) Tort et al. PNAS 2008
 ==== Introduction ====
@@ Line 139: / Line 142: @@
 ☛ Subtract the mean from the rectified signal and plot the acorr again. Verify that the value at lag 1 matches the output of ''corrcoef()'' for lag 1.
 === Application to real data ===
@@ Line 207: / Line 211: @@
 ==== Extracting signal phases ====
+In the previous section, we explored the relationships between different power bands of a signal. Such relationships are a form of //cross-frequency coupling//: what happens in one frequency band depends on what happens in another. Such cross-frequency relationships are not limited to //power//, but can also extend to //phase//, as the following figure (from [[http://www.ncbi.nlm.nih.gov/pubmed/17548233 | Jensen and Colgin, 2007]]) illustrates:
+{{ :analysis:course:jensen_colgin_fig1.png?750 |}}
+The power-to-power coupling (a) we have explored using the correlation matrix of the spectrogram. To reveal cross-frequency coupling involving phase, we need a way to extract the phase of a signal. One way to accomplish this is to use the Hilbert transform:
+<code matlab>
+Fs = 500; dt = 1./Fs;
+tvec = 0:dt:1-dt;
+wnoise = randn(size(tvec));
+sgnl = sin(2*pi*4*tvec) + 0.1*wnoise;
+h = hilbert(sgnl);
+phi = angle(h);
+</code>
+The Hilbert transform of a periodic signal returns a complex quantity that, similar to the Fourier transform, contains both amplitude and phase information. Here, we access the phase component by ''angle(h)'' (the amplitude would be obtained by ''abs(h)''), and plot it with some ''plotyy()'', a nice function that gives two signals their own y-axis:
+<code matlab>
+[ax_h,h1,h2] = plotyy(tvec,sgnl,tvec,phi);
+set(ax_h(2),'YColor',[1 0 0],'YLim',[-pi pi],'YTick',-pi:pi/2:pi,'YTickLabel',{'-p','-p/2','0','p/2','p'}, ...
+   'fontname','symbol','XTick',[]);
+set(h2,'Color',[1 0 0]);
+set(get(ax_h(2),'Ylabel'),'String','phase (radians)')
+box off;
+</code>
+The result:
+{{ :analysis:course:week7_fig6.png?600 |}}
+Obviously, phase is only defined for a signal with sufficiently narrow frequency content. Thus, for real data, the application of the Hilbert transform to extract signal phase should be preceded by a bandpass filter.
+A further caveat when extracting phase is that if the signal gets small relative to the noise, phases still get returned, but they can get unreliable without warning.
+☛ Contaminate the above signal with increasing amounts of noise to verify this.
 ==== Relating signal phases to amplitude ====
+=== Synthetic data ===
+How can we detect relationships between phase and amplitude? It will be helpful to first generate a test signal that displays phase-amplitude coupling:
+<code matlab>
+Fs = 500; dt = 1./Fs;
+tvec = 0:dt:1-dt;
+f1 = 8;
+s1 = sin(2*pi*f1*tvec)+0.05*randn(size(tvec));
+s1_phi = angle(hilbert(s1)); % phase of slow signal
+dphi = exp(-abs(s1_phi)/1.5); % envelope with peaks at phase 0, based on slow phase
+f2 = 80;
+s2 = 0.3.*sin(2*pi*f2*tvec).*dphi; % fast oscillation multiplied by envelope
+s = s1+s2;
+plot(tvec,s); xlim([0 0.5]);
+</code>
+This gives
+{{ :analysis:course:week7_fig7.png?600 |}}
+Notice how the fast (80Hz) oscillation is strongest at the peaks of the slow (8Hz) signal. Given such a signal, we may wonder if we can extract this relationship by computing the correlation coefficient between phase (of the slow signal) and power (of the fast signal). However, as ''corrcoef(s1_phi,dphi)'' will demonstrate, the correlation is in fact zero. This is because the relationship between phase and power is //non-linear//: starting from the most negative phase, power first increases with phase until the phase is zero, and then decreases again in a symmetric fashion. A linear measure such as Pearson's correlation coefficient cannot detect this, as illustrated in the figure in the first section.
+A general-purpose workaround is to divide signal phase into bins, and compute the average power for each bin. We can then ask if power varies as a function of bin, which if true would indicate cross-frequency phase-amplitude coupling.
+To do this, we can define a function ''averageXbyYbin()'' as follows:
+<code matlab>
+function x_avg = averageXbyYbin(x,y,y_edges)
+[~,idx] = histc(y,y_edges); % idx returns which bin each point in y goes into
+x_avg = zeros(size(y_edges));
+for iBin = length(y_edges):-1:1 % for each bin...
+   if sum(idx == iBin) ~= 0 % at least one sample in this bin
+      x_avg(iBin) = nanmean(x(idx == iBin)); % compute average of those x's that go in it
+   end
+end
+</code>
+One quirk of MATLAB's phenomenally useful ''histc()'' function is that it creates N bins out of N edges, when perhaps N-1 bins would be more intuitive. This is because the Nth output bin is for those data points that fall exactly on the Nth edge. Thus, to plot our output correctly, we have to do something like
+<code matlab>
+phi_edges = -pi:pi/8:pi;
+pow_bin = averageXbyYbin(dphi,s1_phi,phi_edges);
+pow_bin(end-1) = pow_bin(end-1)+pow_bin(end); % add counts on last edge to preceding bin
+pow_bin = pow_bin(1:end-1); % trim
+phi_centers = phi_edges(1:end-1)+pi/16; % convert edges to centers
+plot(phi_centers,pow_bin);
+</code>
+The output shows that we have recovered the original phase-amplitude relationship. The ''averageXbyYbin()'' function can be modified to also return standard deviations, or even the full distributions, for each bin to permit statistical significance testing.
+=== Real data ===
+The cross-frequency power correlation matrix above suggests that low-gamma power (~45-65Hz) co-occurs with delta power (3-4Hz). This raises the possibility that low-gamma power is modulated by the phase of the delta rhythm.
+☛ Plot the average low-gamma power as a function of delta phase for the ventral striatal LFP used previously in this module. As with the cross-frequency correlation plot above, restrict your analysis to the risk session only. You can use a skeleton like the following:
+<code matlab>
+%% load and restrict the data
+%% define frequency ranges of interest
+%% design filters for frequency ranges
+%% filter the data (remember to use filtfilt!)
+%% extract delta phase and low gamma power
+%% use averageXbyYbin to plot relationship (ideally with standard deviations)
+</code>
+In fact, the above is unlikely to work "out of the box" because for much of the session, there is very little delta power. Thus, the phases that will be extracted will be essentially random, as we noted above. The code you used to detect sharp wave-ripple complexes in Module 5 can be adapted to exclude epochs without much delta power.

User Tools

Site Tools

Differences

Page Tools