-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathdpac_Fig3b_simulating_asymmetric_phases_bpf.m
executable file
·155 lines (116 loc) · 6.35 KB
/
dpac_Fig3b_simulating_asymmetric_phases_bpf.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
%% dPAC: A method for debiasing phase-amplitude cross-frequency coupling
% Joram van Driel, Roy Cox & Mike X Cohen
% 2014/2015
% --
% This code accompanies the paper titled "Phase clustering bias in
% phase-amplitude cross-frequency coupling and its removal". Below,
% simulations are run that show how non-sinusoidal properties can produce
% phase clustering bias in different coupling measures.
% In addition, band-pass filtering is used to approximate regular time
% series analysis approaches to real data.
% Using the code without following the paper may lead to confusion,
% incorrect data analyses, and misinterpretations of results.
% The authors assume no responsibility for inappropriate or incorrect use
% of this code.
clear all
% cd('some/path/'); % -- change directory if you want to save output and/or plots
%% parameters
srate = 1000; % -- sampling rate
t = 1/srate:1/srate:12; % -- time points: 10000 1ms steps; 12 seconds (10 seconds plus padded 1000 ms on each side for edge artifacts)
ntimepoints = length(t)-2000; % -- get number of timepoints minus padded buffer zones
centTimes = 0:0.2:12; % -- center times for 'oscillatory' peaks; corresponds to 5 Hz
fpow = 30; % -- frequency for power: gamma
nbins = 18; % -- number of bins used for Tort's MI
gausWidth = linspace(0.001,0.08,50); % -- gaussian width of the 'oscillatory' peaks
[pac, dpac, mi, pc, plv] = deal(zeros(2,50)); % -- initialize output matrices
%%
q=0;
for in_anti = 1:2 % -- simulate two scenarios: phase clustering bias as in vs. anti phase with coupling angle
for g = 1:50 % -- loop over 50 widths of the gaussian 'cycles'
gW = gausWidth(g); % -- set the width of the gaussian
% -- create time series
% -- this codes concatenates a series of gaussians that are
% -- then detrended
cmd = 'ts = detrend(';
for i=1:length(centTimes)
cmd = [cmd 'exp(-((t-' num2str(centTimes(i)) ').^2)/(2*' num2str(gW) '^2)) + '];
end
cmd = [cmd '0);'];
eval(cmd)
% -- gamma is a complex sine wave that is phase-modulated by the
% -- above time series
gamma = ((ts+0.5) .* exp( 1i * 2 * pi * fpow * t ));
if in_anti==2
ts = ts([101:end 1:100]); % -- shift the time series with pi (100 ms)
end
thetagamma = real(gamma) + ts;
%% basic filter settings
% -- band-pass filter
thetaband = [3 7];
gammaband = [25 35];
theta_filt_order = round(3*(srate/thetaband(1)));
theta_filterweights = fir1(theta_filt_order,[mean(thetaband)-(thetaband(2)-thetaband(1)) mean(thetaband)+(thetaband(2)-thetaband(1))]/(srate/2)); % -- NOTE: if eeglab is added to path, remove eeglab/plugins/biosig from path
gamma_filt_order = round(3*(srate/gammaband(1)));
gamma_filterweights = fir1(gamma_filt_order,[mean(gammaband)-(gammaband(2)-gammaband(1)) mean(gammaband)+(gammaband(2)-gammaband(1))]/(srate/2)); % -- NOTE: if eeglab is added to path, remove eeglab/plugins/biosig from path
thetafilt = filtfilt(theta_filterweights,1,thetagamma); thetafilt = thetafilt(1001:end-1000);
gammafilt = filtfilt(gamma_filterweights,1,thetagamma); gammafilt = gammafilt(1001:end-1000);
%% PAC, dPAC, MI, PLV (note: dPAC not shown in the paper)
% -- extract phase and power
thetaphase = angle(hilbert(thetafilt));
gammapow = abs(hilbert(gammafilt));
% -- equations to calculate PAC, dPAC and PC (phase clustering vector length)
pac(in_anti,g) = abs(mean(exp(1i*thetaphase) .* gammapow));
debias_term = mean(exp(1i*thetaphase)); % -- this is the phase clustering bias
dpac(in_anti,g) = abs(mean( (exp(1i*thetaphase) - debias_term) .* gammapow)); % -- which is subtracted here
pc(in_anti,g) = abs(mean(exp(1i*thetaphase)));
% -- Tort's Modulation Index
thetaphase_bin = ceil( tiedrank( thetaphase ) / (ntimepoints / nbins) ); % -- bin the theta phase angles into nbins -- NOTE: tiedrank also exists in eeglab toolbox; when added to path, may cause conflict
gammapow_bin = zeros(1,nbins);
for k=1:nbins
gammapow_bin(k) = squeeze(mean(gammapow(thetaphase_bin==k))); % -- compute mean gamma power in each bin
end
gammapow_bin = gammapow_bin ./ sum(gammapow_bin); % -- normalize
mi(in_anti,g) = (log(nbins) + sum(gammapow_bin.*log(gammapow_bin)) ) ./ log(nbins); % -- compute MI
% -- Phase-locking value (Cohen, 2008; Colgin et al 2009)
plv(in_anti,g) = abs(mean(exp(1i*(thetaphase- angle(hilbert(detrend(gammapow)))))));
end
end
%% plot of different CFC measures (Figure 3B)
figure('position',[400 100 400 400])
subplot(221)
plot(gausWidth,pc(1,:),'b','linewidth',1); hold on
plot(gausWidth,pc(2,:),':b','linewidth',2);
set(gca,'xscale','log','xtick',[0.001 0.004 0.016 0.064],'ylim',[0 0.8]) %
title('PC')
box off
% -- small inlay plot to zoom in on small bump of PC increase for PC being
% -- anti-phase with CFC angle (not visible because of strong PC increase
% -- for low-g and PC-CFC in-phase
ax1 = gca;
ax1_pos = get(ax1,'position');
ax2 = axes('Position',[0.27 0.78 0.15 0.1],...
'XAxisLocation','bottom',...
'YAxisLocation','left',...
'Color','none');
plot(gausWidth,pc(1,:),'parent',ax2,'color','b','linewidth',1); hold on
set(gca,'xscale','log','xtick',[0.001 0.064],'ylim',[0 0.006])
box off
subplot(222)
plot(gausWidth,pac(1,:),'r','linewidth',1); hold on
plot(gausWidth,pac(2,:),':r','linewidth',2);
set(gca,'xscale','log','xtick',[0.001 0.004 0.016 0.064])
title('(d)PAC')
box off
legend('anti-phase','in-phase')
subplot(223)
plot(gausWidth,mi(1,:),'r','linewidth',1); hold on
plot(gausWidth,mi(2,:),':r','linewidth',2);
set(gca,'xscale','log','xtick',[0.001 0.004 0.016 0.064],'ylim',[0 0.1])
title('MI')
box off
subplot(224)
plot(gausWidth,plv(1,:),'r','linewidth',1); hold on
plot(gausWidth,plv(2,:),':r','linewidth',2);
set(gca,'xscale','log','xtick',[0.001 0.004 0.016 0.064])
title('PLV')
box off