Robust circular cluster-based statistics for respiration-brain coupling.

1. %% Data preparation script: respiratory phase extraction, surrogate generation & hitrate binning
2. %
3. % This script extracts respiratory phase from raw respiration data using
4. % two-point interpolation, generates surrogate phases via the Iterative
5. % Amplitude-Adjusted Fourier Transform (IAAFT), and bins behavioral data
6. % (e.g., hit rates) into empirical and surrogate respiratory phase bins.
7. % A detailed description of the preparation pipeline can be found in
8. % [ref].
9. %
10. % Input:
11. % - raw respiration time series vector
12. % - events structure with stimulus onset and outcome field for each subject
13. % .onset: contains timestamp of stimulus presentation for each trial (in samples from start of recording)
14. % .outcome: contains outcome value for each trial (here: hit vs miss)
15. %
16. % Output (here: hit rate):
17. % binhr.mat = mean empirical outcome for each respiratory phase bin, for each subject
18. % (subjects x phase bins matrix)
19. % sbinhr.mat = mean surrogate outcome for each respiratory phase bin and surrogate, for each subject
20. % (subjects x phase bins x surrogate iterations matrix)
21. %
22. % Required Helper Functions:
23. % - two_point_interp: extracts phase of a signal by two-point interpolation
24. % - generate_surrogate_iaaft: creates a phase-scrambled surrogate time series
25. % of a signal
26. % - PLV: computes phase-locking value of two phase time series
27. %
28. %
29. % Copyright (C) 2025, Daniel Kluger & Teresa Berther, University of Münster, Germany
31. clearvars
32. clc
33. close all
35. addpath(genpath('~/respmethods/matlab/'));
36. datapath = '~/respmethods/_exampledata/'; % contains raw respiration traces and events table
39. % subject ids
40. ids = {'001' '002' '003' '004' '005' '006' '007' '009' '010' '011' '012' '013' '014' '015' '016' '019' '020' '022' '023'};
42. % config
43. orgfs = 600; % original sampling frequency of raw data
44. fs = 100; % final sampling frequency after downsampling
46. niter = 5000; % number of surrogate iterations
47. plvcrit = 0.1; % threshold for IAAFT surrogate generation, max PLV of surrogate with data
48. nbin = 60; % number of phase bins
49. phw = (2*pi/10)/2; % width of each phase bin, divided in half
50. pb = linspace(-pi,pi,nbin+1); % vector containing centre of each phase bin (= phase bin vector)
51. pb(end) = [];
53. load(fullfile(datapath, 'events.mat')); % load stimulus onsets ("onsets") and outcome (here: "corr" for hit vs miss trials)
55. %% Compute true and surrogate respiration phase values at stimulus onset
56. for isub = 1:numel(ids)
58. disp(['Computing phase and surrogates for subject #' num2str(isub) '/' num2str(length(ids))]);
60. pattern = sprintf('*_%s.mat', ids{isub});
61. respfiles = dir(fullfile(datapath, pattern));
62. load([datapath respfiles.name]);
64. x = downsample(zscore(data(1,:)), orgfs/fs); % downsample & zscore raw respiration signal
65. x = movmean(x, 0.4*fs); % some smoothing
67. % extract phase vector of empirical data using two-point interpolation
68. [pv, first, last] = two_point_interp(x); % first/last indices saved for NaN exclusion later
70. %%% generate surrogate phase vectors using IAAFT
71. allsresp = nan(niter,length(x));
72. k = 1;
73. while k <= niter
75. disp(['Computing surrogate resp trace #' num2str(k) '/' num2str(niter)]);
77. sresp = generate_surrogate_iaaft(x,'verbose',false); % single mock solution for respiration time series
78. [spv, ftmp, ltmp] = two_point_interp(sresp); % extract surrogate phase via two-point interpolation
79. spvtmp = spv(max(first, ftmp):min(last, ltmp)); % temporarily get rid of NaNs in phase vectors as PLV function cannot handle NaNs
80. pvtmp = pv(max(first, ftmp):min(last, ltmp));
81. plv = PLV(pvtmp', spvtmp); % compute phase-locking value
82. if plv < plvcrit % only keep solution if plv is smaller than current constraint
83. allsresp(k,:) = spv; % save full surrogate phase vector
84. k = k+1;
85. else
86. continue % otherwise start over
87. end
88. end
89. save(fullfile(datapath, ['surrogates_iaaft_' num2str(isub) '.mat']),'allsresp', '-v7.3');
91. % get stimulus onsets in the respiration time frame
92. % (= time from recording start in ms)
93. disp(['Computing empirical and surrogate onsets for subject #' num2str(isub) '/' num2str(length(ids))]);
95. onsets = events{isub}.onsets; % the onsets field contains the sample time stamps for the stimuli presentation
96. onsets = round(onsets/(orgfs/fs)); % downsample the sample time stamps to match new fs
97. empphase = pv(onsets); % empirical phase at onset of stimulus for each trial
98. surrphases = [];
99. for k = 1:niter
100. surrphases(:,k) = allsresp(k, onsets); % surrogate phase at onset of stimulus for each trial
101. end
103. % save
104. save(fullfile(datapath, ['empiricalphase_' num2str(isub) '.mat']),'empphase');
105. save(fullfile(datapath, ['surrogatephases_' num2str(isub) '.mat']),'surrphases');
107. end
110. %% Compute outcome (hit rates) ~ respiration phase
111. % example outcome here is hit rate, but can be any outcome variable
112. binhr = nan(numel(ids), nbin);
113. sbinhr = nan(numel(ids), nbin, niter);
115. for isub = 1:numel(ids)
117. disp(['Computing empirical and surrogate binned hit rates for subject #' num2str(isub) '/' num2str(length(ids))]);
119. % empirical binHR
120. load(fullfile(datapath, ['empiricalphase_' num2str(isub) '.mat']),'empphase');
121. for ibin = 1:nbin
122. phsel = find((empphase>pb(ibin)-phw) & (empphase<pb(ibin)+phw) | ... % find all trials whose phase falls within the current phase bin
123. (empphase-2*pi>pb(ibin)-phw) & (empphase-2*pi<pb(ibin)+phw)| ...
124. (empphase+2*pi>pb(ibin)-phw) & (empphase+2*pi<pb(ibin)+phw));
125. binhr(isub,ibin) = sum(events{isub}.outcome(phsel))/numel(phsel); % based on the outcome for these trials, compute hit rate for this phase bin
126. end
128. % surrogate binHR
129. load(fullfile(datapath, ['surrogatephases_' num2str(isub) '.mat']),'surrphases');
130. for k = 1:niter
131. sphase = surrphases(:, k);
132. for ibin = 1:nbin
133. phsel = find((sphase>pb(ibin)-phw) & (sphase<pb(ibin)+phw) | ... % again find trials for each phase bin...
134. (sphase-2*pi>pb(ibin)-phw) & (sphase-2*pi<pb(ibin)+phw)| ...
135. (sphase+2*pi>pb(ibin)-phw) & (sphase+2*pi<pb(ibin)+phw));
136. sbinhr(isub,ibin,k) = sum(events{isub}.outcome(phsel))/numel(phsel); % ...and compute bin HR based on the outcome
137. end
138. end
139. end
141. % save mean hit rates ~ respiration phase bins for empirical & surrogate data
142. % this is the data we run the circular clustering on
143. save(fullfile(datapath, 'binhr.mat'), 'binhr'); % subjects x nbins array of empirical phase-binned outcome values
144. save(fullfile(datapath, 'surrbinhr.mat'), 'sbinhr'); % subjects x nbins x niter array of surrogate phase-binned values
145. save(fullfile(datapath, 'phasebinvect.mat'), 'pb'); % also save the vector with the phase bin centres

Tutorial_DataPrep.m at commit f2bd9a0, no license · at the source