event_detection/example_script.m at main · SPiN-Lab/event_detection · GitHub

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
clear all
close all
clc

load ts

% zscore time series
z = zscore(ts);

% number of time points/nodes
[t,n] = size(z);

% upper triangle indices (node pairs = edges)
[u,v] = find(triu(ones(n),1));

% edge time series
ets = z(:,u).*z(:,v);

% calculate rss
rss = sum(ets.^2,2).^0.5;

% repeat with randomized time series
numrand = 100;

% initialize array for null rss
rssr = zeros(t,numrand);

% perform numrand randomizations
for irand = 1:numrand

    % create circularly shifted time series
    zr = z;
    for i = 1:n
        zr(:,i) = circshift(zr(:,i),randi(t));
    end

    % edge time series with circshift data
    etsr = zr(:,u).*zr(:,v);

    % calcuate rss
    rssr(:,irand) = sum(etsr.^2,2).^0.5;

end

% calculate p-value
p = zeros(t,1);
for i = 1:t
    p(i) = mean(rssr(:) >= rss(i));
end

% apply statistical cutoff
pcrit = 0.001;

% find frames that pass statistical test
idx = find(p < pcrit);

% identify contiguous segments of frames that pass statistical test
dff = idx' - (1:length(idx));
unq = unique(dff);
nevents = length(unq);

% find the peak rss within each segment
idxpeak = zeros(nevents,1);
for ievent = 1:nevents
    idxevent = idx(dff == unq(ievent));
    rssevent = rss(idxevent);
    [~,idxmax] = max(rssevent);
    idxpeak(ievent) = idxevent(idxmax);
end

% get activity at peak
tspeaks = z(idxpeak,:);

% get co-fluctuation at peak
etspeaks = tspeaks(:,u).*tspeaks(:,v);

%% plot rss time series, null, and significant peaks
figure('position',[200,200,1400,200])
ph = plot(1:t,rssr,'color',ones(1,3)*0.65);
hold on;
qh = plot(idxpeak,rss(idxpeak),'r*',1:t,rss,'k','linewidth',2);
xlim([1,t])
xlabel('frame'); ylabel('rss');
legend([ph(1); qh],'null','significant','orig');

%% calculate mean co-fluctuation (edge time series) across all peaks
mu = nanmean(etspeaks,1);

% represent in matrix form
mat = zeros(n);
mat(triu(ones(n),1) > 0) = mu;
mat = mat + mat';

% load brain systems from Gordon et al
load hcp333
[~,idxsort] = sort(lab);

% draw matrix of co-fluctuation magnitude
figure('position',[200,200,600,600]);
axes('outerposition',[0.15,0.15,0.85,0.85]);
imagesc(mat(idxsort,idxsort),[-3,3]);
axis off

% add lines between systems
hold on;
idx = find(diff(lab(idxsort)));
for j = 1:length(idx)
    plot([0.5,n + 0.5],(idx(j) + 0.5)*ones(1,2),'k')
    plot((idx(j) + 0.5)*ones(1,2),[0.5,n + 0.5],'k')
end

% add system names
for i = 1:max(lab)
    x = mean(find(lab(idxsort) == i));
    text(-0.01*n,x,net{i},'horizontalalignment','right')
    text(x,1.01*n,net{i},'horizontalalignment','right','rotation',90)
end