OSCR

Theta oscillations are an organizational unit of odor processing in the olfactory bulb.

Code ↔ Paper

16 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 16 matches
  1. [1] § MATERIALS AND METHODS › Cross-frequency coupling of theta and gamma ↔ code_fig6.m, lines 426–474 · score 0.95 · Hilbert transform, gamma bursts, 30–60 Hz, gamma peaks, filtered signal, CircStat
  2. [2] § MATERIALS AND METHODS › Cross-frequency coupling of theta and gamma ↔ dataScript.m, lines 498–533 · score 0.94 · low gamma filtered, 30–60 Hz, amplitude envelope, gamma envelope, low gamma frequency, filtered signal
  3. [3] § MATERIALS AND METHODS › Intertrial phase coherence ↔ code_suppFigs.m, lines 543–680 · score 0.83 · respiratory phase bins, circularly shifted, passive inhalations, respiratory epochs, intentional sniffs, Fisher
  4. [4] § RESULTS › Theta oscillations organize the timing of gamma bursts in the olfactory bulb ↔ code_fig6.m, lines 426–474 · score 0.80 · gamma burst, 30–60 Hz, gamma peaks, theta phase, low gamma, sniff onset
  5. [5] § RESULTS › Theta oscillations organize the timing of gamma bursts in the olfactory bulb ↔ code_fig6.m, lines 301–424 · score 0.77 · aligned spectrogram, theta troughs, Theta peak, theta filtered, Amplitude spectrograms, theta frequency
  6. [6] § MATERIALS AND METHODS › Intertrial phase coherence ↔ code_fig4.m, lines 585–639 · score 0.77 · respiratory phase bins, circularly shifted, phase epochs, ITPC, matched, permuted
  7. [7] § MATERIALS AND METHODS › Statistical analysis ↔ code_fig4.m, lines 201–304 · score 0.71 · sign flipped, cluster mass, clusters formed, intentional sniffs, passive, threshold
  8. [8] § RESULTS › High-precision recordings from the human olfactory bulb ↔ code_fig3.m, lines 95–186 · score 0.67 · low gamma band, Low gamma odor, odor detection, High gamma, 150 Hz, 60 Hz
  9. [9] § MATERIALS AND METHODS › Amplitude spectrograms ↔ code_suppFigs.m, lines 403–541 · score 0.66 · temporally smoothed, longest, shortest, S4C, bandwidth, SD
  10. [10] § RESULTS › Intentional sniffs enhance and align theta oscillations in the olfactory bulb ↔ code_fig4.m, lines 201–304 · score 0.57 · cluster corrected permutation, passive inhales, intentional sniffs, theta amplitude, airflow, bar
  11. [11] § RESULTS › Sniffing behavior covaries with olfactory bulb theta oscillations ↔ code_suppFigs.m, lines 403–541 · score 0.55 · sniff duration, passive inhales, intentional sniffs, S4C, subset, overlapping
  12. [12] § RESULTS › Intentional sniffs enhance and align theta oscillations in the olfactory bulb ↔ code_fig4.m, lines 776–805 · score 0.55 · maximum theta ITPC, passive inhales, intentional sniff, resampling, phase, Figure 4
  13. [13] § RESULTS › High-precision recordings from the human olfactory bulb ↔ code_fig2_OdorImageCue.m, lines 329–393 · score 0.55 · image onset, odor cues, FDR corrected, baseline, Figure 2, threshold
  14. [14] § RESULTS › High-precision recordings from the human olfactory bulb ↔ code_fig3.m, lines 274–372 · score 0.53 · Power spectra, Gamma power, cue task, low gamma, contact, Figure 3
  15. [15] § RESULTS › High-precision recordings from the human olfactory bulb ↔ code_fig2_OdorImageCue.m, lines 329–393 · score 0.53 · image cues, Odor cues, FDR corrected, baseline, onset, Figure 2
  16. [16] § RESULTS › Intentional sniffs enhance and align theta oscillations in the olfactory bulb ↔ code_suppFigs.m, lines 682–719 · score 0.51 · Resampled distributions, passive inhales, intentional sniffs, CI, bootstrap, subset

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

MATLAB · 1,002 lines · 45 KB · CC-BY-4.0 · 4 matches

  1. %% Load data
  2. load data_supp
  4. %% Fig S2A - Odor vs. Image spectrograms
  5. % CUE TASK ODOR
  6. % Calculate spectrogram - concatenated 6 subs
  7. tic
  8. ts = data_supp.lfpData_cueSniffOdor_combinedSubs;
  9. sniffonset_conds = { data_supp.sniffonsets_cueOdorSniff_combinedSubs };
  10. para = [];
  11. para.freq = logspace( log10( 1), log10( 200), 100); % center frequencies for filtering
  12. para.bandwidth = logspace( log10( 2), log10( 30), 100); % bandwidth of each center frequency
  13. para.srate = srate; % sampling rate, in Hz
  14. para.events = sniffonset_conds; % in samples
  15. para.epoch = [-1, 5]; % in seconds, relative to events
  16. para.baseline = [-0.6, -0.1]; % in seconds, relative to events
  17. para.smooth_win = 0.01; % temporal smoothing, in secondsj
  18. para.nbpermuts = 10000; % number of permutations,
  19. para.measure = 'amplitude'; % or 'power'
  20. para.baseline_condition = []; % which condition to use as baseline (empty for just before each of multiple conditions)
  21. para.trl = 'mean'; % 'mean' o r'median' value across trials
  22. [stats_cueSniffOdor, epoch_time_cueSniffOdor, ev_epochs_cueSniffOdor, ~, sd_cueSniffOdor] = PermutAmpChange_Commonbaseline(ts, para);
  24. % IMAGE ONSETS
  25. ts = data_supp.lfp_visualCue_combinedSubs;
  26. sniffonset_conds = { data_supp.imageonsets_combinedSubs };
  27. para = [];
  28. para.freq = logspace( log10( 1), log10( 200), 100); % center frequencies for filtering
  29. para.bandwidth = logspace( log10( 2), log10( 30), 100); % bandwidth of each center frequency
  30. para.srate = srate; % sampling rate, in Hz
  31. para.events = sniffonset_conds; % in samples
  32. para.epoch = [-1, 5]; % in seconds, relative to events
  33. para.baseline = [-0.6, -0.1]; % in seconds, relative to events
  34. para.smooth_win = 0.01; % temporal smoothing, in secondsj
  35. para.nbpermuts = 10000; % number of permutations,
  36. para.measure = 'amplitude'; % or 'power'
  37. para.baseline_condition = []; % which condition to use as baseline (empty for just before each of multiple conditions)
  38. para.trl = 'mean'; % 'mean' o r'median' value across trials
  39. [stats_imageOnsets, epoch_time_imageOnsets, ev_epochs_imageOnsets, ~, sd_imageOnsets] = PermutAmpChange_Commonbaseline(ts, para);
  41. % CALCULATE DIFFERENCE BETWEEN CUE TASK SNIFF ODOR AND CUE TASK VISUAL IMAGE
  42. mat = cat(3,ev_epochs_cueSniffOdor{1},ev_epochs_imageOnsets{1});
  43. mat = permute(mat, [3 1 2]);
  44. g = [ones(length(data_supp.sniffonsets_cueOdorSniff_combinedSubs),1) ; 2*ones(length(data_supp.imageonsets_combinedSubs),1)];
  45. AmpChange_CueTasks_OdorVsImage = G_PowChange_Clust(mat, g, epoch_time_imageOnsets, para.freq, 'baseline', para.baseline,'db_transf', 'no','common_bs','no');
  47. % PLOT DIFFERENCE BETWEEN CUE TASK SNIFF ODOR AND CUE TASK VISUAL IMAGE
  48. % Respiratory overlays
  49. resp_epochs_cueSniffOdor_combinedSubs = squeeze(zscore(G_EpochTS(data_supp.resp_cueSniffOdor_combinedSubs, data_supp.sniffonsets_cueOdorSniff_combinedSubs, [-1,5]*srate)));
  50. resp_meanEpochs_cueSniffOdor_combinedSubs = squeeze(mean(resp_epochs_cueSniffOdor_combinedSubs,2));
  51. resp_epochs_imageOnsets_combinedSubs = squeeze(zscore(G_EpochTS(data_supp.resp_visualCue_combinedSubs, data_supp.imageonsets_combinedSubs, [-1,5]*srate)));
  52. resp_meanEpochs_imageOnsets_combinedSubs = squeeze(mean(resp_epochs_imageOnsets_combinedSubs,2));
  54. % Calculate significance threshld using FDR correction
  55. %Odor Cue Task - Sniff Odor
  56. P_cueSniffOdor.P = G_z2p(stats_cueSniffOdor.Z{1});
  57. FDR_cueSniffOdor_Odor.FDR = fdr(P_cueSniffOdor.P(:));
  58. thresh_cueSniffOdor_Odor.thresh = G_p2z(FDR_cueSniffOdor_Odor.FDR);
  59. %Visual Cue Task - Image Onset
  60. P_visualCue_image.P = G_z2p(stats_imageOnsets.Z{1});
  61. FDR_visualCue_image.FDR = fdr(P_visualCue_image.P(:));
  62. thresh_visualCue_image.thresh = G_p2z(FDR_visualCue_image.FDR);
  63. %amplitude difference
  64. AmpChange_P= G_z2p(AmpChange_CueTasks_OdorVsImage.statsmap);
  65. AmpChange_FDR = fdr(AmpChange_P(:));
  66. AmpChange_thresh = G_p2z( AmpChange_FDR);
  68. % Plot spectrogram
  69. Time = -1 : 1/srate : 5;
  70. subplot(3,1,1)
  71. % subplot(1,3,1)
  72. G_Imagesc(stats_cueSniffOdor.Z{1},'yaxis',para.freq,'overlay',resp_meanEpochs_cueSniffOdor_combinedSubs,...
  73. 'xaxis', Time, 'xtick', [-0.3,0:0.5:5], 'threshold',thresh_cueSniffOdor_Odor.thresh,...
  74. 'ytick', [1,3,8,15,30,50,100,200],'Fontsize',14,'clim',[-10,10],'ctick', [-10,0,10],...
  75. 'ytick_precision',0,'xlabel','Time (s)','ylabel','Frequency (Hz)');
  76. xlim([0.69*srate,3.51*srate]);
  77. title('Odor');
  78. subplot(3,1,2)
  79. % subplot(1,3,2)
  80. G_Imagesc(stats_imageOnsets.Z{1},'yaxis',para.freq,'overlay',resp_meanEpochs_imageOnsets_combinedSubs,...
  81. 'xaxis', Time, 'xtick', [-0.3,0:0.5:5], 'threshold',thresh_visualCue_image.thresh,...
  82. 'ytick', [1,3,8,15,30,50,100,200],'Fontsize',14,'clim',[-10,10],'ctick', [-10,0,10],...
  83. 'ytick_precision',0,'xlabel','Time (s)','ylabel','Frequency (Hz)');
  84. xlim([0.69*srate,3.51*srate]);
  85. title('Image');
  86. subplot(3,1,3)
  87. % subplot(1,3,3)
  88. h=G_Imagesc(AmpChange_CueTasks_OdorVsImage.statsmap,'yaxis',para.freq,...
  89. 'xaxis', Time, 'xtick', [-0.3,0:0.5:5], 'threshold',AmpChange_thresh,...
  90. 'ytick', [1,3,8,15,30,50,100,200],'Fontsize',14,...
  91. 'clim',[-5,5],'ctick', [-5,0,5],...
  92. 'ytick_precision',0,'xlabel','Time (s)','ylabel','Frequency (Hz)');
  93. xlim([0.69*srate,3.51*srate]);
  94. title('Odor minus Image');
  95. % for k = 1:length(AmpChange_CueTasks_OdorVsImage.clust{1}.both.info.PixelIdxList)
  96. % mask = zeros(AmpChange_CueTasks_OdorVsImage.clust{1}.both.info.ImageSize);
  97. % if AmpChange_CueTasks_OdorVsImage.clust{1}.both.statssum_pval(k) < 0.05
  98. % mask(AmpChange_CueTasks_OdorVsImage.clust{1}.both.info.PixelIdxList{k}) = 1;
  99. % G_Outline_Pixels(h,mask,'linecolor','m','linewidth',1);
  100. % end
  101. % end
  102. % set(gcf,'pos',[367 571 1722 411]); % For 3 columns x 1 row
  103. set(gcf,'pos',[1000 281 507 957]); % For 1 column x 3 rows
  104. xlim([0.69*srate,3.51*srate]);
  105. ax = gca; ax.LineWidth = 1; ax.FontSize = 16;
  107. %% Fig S2A - Image onsets no sniff
  108. ts = data_supp.lfp_visualCue_combinedSubs;
  109. sniffonset_conds = { data_supp.imageonsets_noSniff_combinedSubs };
  110. para = [];
  111. para.freq = logspace( log10( 1), log10( 200), 100); % center frequencies for filtering
  112. para.bandwidth = logspace( log10( 2), log10( 30), 100); % bandwidth of each center frequency
  113. para.srate = srate; % sampling rate, in Hz
  114. para.events = sniffonset_conds; % in samples
  115. para.epoch = [-1, 5]; % in seconds, relative to events
  116. para.baseline = [-0.6, -0.1]; % in seconds, relative to events
  117. para.smooth_win = 0.01; % temporal smoothing, in secondsj
  118. para.nbpermuts = 10000; % number of permutations,
  119. para.measure = 'amplitude'; % or 'power'
  120. para.baseline_condition = []; % which condition to use as baseline (empty for just before each of multiple conditions)
  121. para.trl = 'mean'; % 'mean' o r'median' value across trials
  122. [stats_imageOnsetsNoSniff, epoch_time_imageOnsetsNoSniff, ev_epochs_imageOnsetsNoSniff, ~, sd_imageOnsetsNoSniff] = PermutAmpChange_Commonbaseline(ts, para);
  124. % Respiratory overlay
  125. resp_epochs_imageOnsetsNoSniff_combinedSubs = squeeze(zscore(G_EpochTS(data_supp.resp_visualCue_combinedSubs, data_supp.imageonsets_noSniff_combinedSubs, [-1,5]*srate)));
  126. resp_meanEpochs_imageOnsetsNoSniff_combinedSubs = squeeze(mean(resp_epochs_imageOnsetsNoSniff_combinedSubs,2));
  128. % Significance threshold
  129. P_imageOnsetsNoSniff.P = G_z2p(stats_imageOnsetsNoSniff.Z{1});
  130. FDR_imageOnsetsNoSniff.FDR = fdr(P_imageOnsetsNoSniff.P(:));
  131. thresh_imageOnsetsNoSniff.thresh = G_p2z(FDR_imageOnsetsNoSniff.FDR);
  133. % Plot spectrogram
  134. G_Imagesc(stats_imageOnsetsNoSniff.Z{1},'yaxis',para.freq,'overlay',resp_meanEpochs_imageOnsetsNoSniff_combinedSubs,...
  135. 'xaxis', Time, 'xtick', [-0.3,0:0.5:5], 'threshold',thresh_imageOnsetsNoSniff.thresh,...
  136. 'ytick', [1,3,8,15,30,50,100,200],'Fontsize',14,'clim',[-10,10],'ctick', [-10,0,10],...
  137. 'ytick_precision',0,'xlabel','Time (s)','ylabel','Frequency (Hz)');
  139. %% Fig S3 - Difference spectrogram for odor vs. air on detection task
  141. % Calculate Subject-wiser Spectrograms at Gamma Frequencies
  142. for sub = 1:num_subs
  143. ts = data_supp.lfp_detect{sub};
  144. sniffonset_conds = { data_supp.sniffonsets_odor{sub}, data_supp.sniffonsets_air{sub} };
  145. para = [];
  146. % para.freq = logspace( log10( 1), log10( 200), 100); % center frequencies for filtering
  147. % para.bandwidth = logspace( log10( 2), log10( 30), 100); % bandwidth of each center frequency
  148. para.freq = linspace( 35, 175, 40); % center frequencies for filtering
  149. para.bandwidth = linspace( 10, 35, length(para.freq)); % bandwidth of each center frequency
  150. para.srate = srate; % sampling rate, in Hz
  151. para.events = sniffonset_conds; % in samples
  152. para.epoch = [-1, 5]; % in seconds, relative to events
  153. para.baseline = [-0.6, -0.1]; % in seconds, relative to events
  154. para.smooth_win = 0.01; % temporal smoothing, in secondsj
  155. para.nbpermuts = 10000; % number of permutations,
  156. para.measure = 'amplitude'; % or 'power'
  157. para.baseline_condition = 1; % which condition to use as baseline (empty for just before each of multiple conditions)
  158. para.trl = 'mean'; % 'mean' o r'median' value across trials
  159. [stats{sub}, epoch_time{sub}, ev_epochs{sub}, ~, sd{sub}] = PermutAmpChange_Commonbaseline(ts, para);
  160. end
  162. z_odor = []; z_air = [];
  163. for sub = 1:num_subs
  164. z_odor(sub,:,:) = stats{sub}.Z{1};
  165. z_air(sub,:,:) = stats{sub}.Z{2};
  166. end
  168. % Average spectrogram output across subjects
  169. z_mean_odor = squeeze(mean(z_odor,1));
  170. z_mean_air = squeeze(mean(z_air,1));
  172. % Calculate amplitude difference between Passive inhale vs. Cued sniff combined subjects
  173. % Combine ev_epochs across subjects
  174. ev_epochs_6normosmics_odor = cat( 3,ev_epochs{1}{1},ev_epochs{2}{1},ev_epochs{3}{1},ev_epochs{4}{1},ev_epochs{5}{1},ev_epochs{6}{1} );
  175. ev_epochs_6normosmics_air = cat( 3,ev_epochs{1}{2},ev_epochs{2}{2},ev_epochs{3}{2},ev_epochs{4}{2},ev_epochs{5}{2},ev_epochs{6}{2} );
  177. mat = cat(3,ev_epochs_6normosmics_odor,ev_epochs_6normosmics_air);
  178. mat = permute(mat, [3 1 2]);
  179. g = [ones(size(ev_epochs_6normosmics_odor,3),1) ; 2*ones(size(ev_epochs_6normosmics_air,3),1)];
  181. AmpChange_odorVair_comBase = G_PowChange_Clust(mat, g, epoch_time{sub}, ...
  182. para.freq, 'baseline', para.baseline,'db_transf', 'no','common_bs','yes');
  184. % Significance threshold
  185. AmpChange_P = []; AmpChange_FDR=[]; AmpChange_thresh=[];
  186. AmpChange_P = G_z2p(AmpChange_odorVair_comBase.statsmap);
  187. AmpChange_FDR = fdr(AmpChange_P(:));
  188. AmpChange_thresh = G_p2z( AmpChange_FDR)
  190. % Plot difference spectrogram
  191. h=G_Imagesc(AmpChange_odorVair_comBase.statsmap,'yaxis',para.freq,...
  192. 'xaxis', Time, 'xtick', [-0.3,0:0.5:5],'threshold',AmpChange_thresh,...
  193. 'ytick', 35:10:175,'Fontsize',14,...
  194. 'clim',[-4,4],'ctick', [-4,-2,0,2,4],...
  195. 'ytick_precision',0);
  196. for k = 1:length(AmpChange_odorVair_comBase.clust{1}.both.info.PixelIdxList)
  197. mask = zeros(AmpChange_odorVair_comBase.clust{1}.both.info.ImageSize);
  198. if AmpChange_odorVair_comBase.clust{1}.both.statssum_pval(k) < 0.05
  199. mask(AmpChange_odorVair_comBase.clust{1}.both.info.PixelIdxList{k}) = 1;
  200. G_Outline_Pixels(h,mask,'linecolor','k','LineWidth',1.5);
  201. end
  202. end
  203. xlim([0.69*srate,3.51*srate]);
  204. title('Difference');
  205. set(gcf,'pos',[334 553 1407 299]);
  207. %% Fig S4 - controlling for sniff features
  208. num_subs=6;
  209. numTrials_passive = cellfun(@length,data_supp.exhale_onsets_passive_bySub);
  210. numTrials_cuedSniff = cellfun(@length,data_supp.exhale_onsets_cuedSniff_bySub);
  211. inhaleDur_passive=[];inhaleDur_cuedSniff=[]; inhalePk_CuedSniff_bySub=[];
  212. time2inhalePk_ind_CuedSniff_bySub=[]; inhalePk_Passive_bySub=[];
  213. time2inhalePk_ind_Passive_bySub=[]; time2inhalePk_ms_Passive=[]; time2inhalePk_ms_CuedSniff=[];
  214. for sub = 1:num_subs
  215. inhaleDur_passive{sub} = (data_supp.exhale_onsets_passive_bySub{sub} - data_supp.sniffonsets_Passive{sub}) / srate;
  216. inhaleDur_cuedSniff{sub} = (data_supp.exhale_onsets_cuedSniff_bySub{sub} - data_supp.sniffonsets_CuedSniff{sub}) / srate;
  217. for i = 1:numTrials_cuedSniff(sub)
  218. [inhalePk_CuedSniff_bySub{sub}(i,1), time2inhalePk_ind_CuedSniff_bySub{sub}(i,1)] = max(resp_epochs_CuedSniff_bySub{sub}(500:1000,i));
  219. end
  220. for i = 1:numTrials_passive(sub)
  221. [inhalePk_Passive_bySub{sub}(i,1), time2inhalePk_ind_Passive_bySub{sub}(i,1)] = max(resp_epochs_Passive_bySub{sub}(500:1000,i));
  222. end
  223. time2inhalePk_ms_Passive{sub} = time2inhalePk_ind_Passive_bySub{sub} / srate * 1000;
  224. time2inhalePk_ms_CuedSniff{sub} = time2inhalePk_ind_CuedSniff_bySub{sub} / srate * 1000;
  225. end
  227. %% Fig S4B - histograms of sniff features
  228. % Check how variables look with histograms for each subject
  229. % Passive
  230. tiledlayout(3,2)
  231. for sub = 1:num_subs
  232. nexttile
  233. histogram(inhaleDur_passive{sub})
  234. % histogram(inhalePk_Passive_bySub{sub})
  235. % histogram(time2inhalePk_ms_Passive{sub})
  236. sgtitle('passive');
  237. end
  238. figure
  239. for sub = 1:num_subs
  240. nexttile
  241. histogram(inhaleDur_cuedSniff{sub})
  242. % histogram(inhalePk_CuedSniff_bySub{sub})
  243. % histogram(time2inhalePk_ms_CuedSniff{sub})
  244. sgtitle('intentional');
  245. end
  247. sub_titles = {'P1','P2','P3','P4','P5','P6'};
  248. tiledlayout(3,2)
  249. for sub = 1:num_subs
  250. nexttile
  251. histogram(inhaleDur_passive{sub},'FaceColor',[1 0.6 0])
  252. hold on
  253. histogram(inhaleDur_cuedSniff{sub},'FaceColor',[0 0.6 1])
  254. if sub == 1
  255. ylabel('Trial count');
  256. legend({'Passive','Intentional'},'Location','northwest');
  257. end
  258. if sub == 5
  259. xlabel('Inhale duration (s)');
  260. end
  261. ax=gca; ax.FontSize = 14; ax.TickDir = 'out'; ax.LineWidth = 1; box off
  262. title(sub_titles{sub});
  264. end
  265. set(gcf,'pos',[300 300 859 702]);
  267. %% Fig S4A - summary statistics of sniff kinematics
  268. % Data structure for mixed effects model
  269. nSub = 6;
  271. DV1_all = [];
  272. DV2_all = [];
  273. DV3_all = [];
  274. Condition = [];
  275. Subject = [];
  277. for s = 1:nSub
  279. % ===== Condition 1 =====
  280. dv1 = inhaleDur_passive{s}(:);
  281. dv2 = inhalePk_Passive_bySub{s}(:);
  282. dv3 = time2inhalePk_ms_Passive{s}(:);
  284. % Sanity check (important!)
  285. assert(numel(dv1)==numel(dv2) && numel(dv1)==numel(dv3), ...
  286. 'Mismatch in trial counts for subject %d, condition C1', s);
  288. n = numel(dv1);
  290. DV1_all = [DV1_all; dv1];
  291. DV2_all = [DV2_all; dv2];
  292. DV3_all = [DV3_all; dv3];
  294. Condition = [Condition; repmat("C1", n, 1)];
  295. Subject = [Subject; repmat(s, n, 1)];
  297. % ===== Condition 2 =====
  298. dv1 = inhaleDur_cuedSniff{s}(:);
  299. dv2 = inhalePk_CuedSniff_bySub{s}(:);
  300. dv3 = time2inhalePk_ms_CuedSniff{s}(:);
  302. assert(numel(dv1)==numel(dv2) && numel(dv1)==numel(dv3), ...
  303. 'Mismatch in trial counts for subject %d, condition C2', s);
  305. n = numel(dv1);
  307. DV1_all = [DV1_all; dv1];
  308. DV2_all = [DV2_all; dv2];
  309. DV3_all = [DV3_all; dv3];
  311. Condition = [Condition; repmat("C2", n, 1)];
  312. Subject = [Subject; repmat(s, n, 1)];
  313. end
  315. % Build the table
  316. tbl = table;
  317. tbl.Subject = categorical(Subject);
  318. tbl.Condition = categorical(Condition);
  319. tbl.duration = DV1_all;
  320. tbl.peak = DV2_all;
  321. tbl.timeToPeak = DV3_all;
  322. tbl.Condition = reordercats(tbl.Condition,{'C1','C2'});
  324. % Sanity checks
  325. height(tbl)
  326. varfun(@numel,tbl,'InputVariables','duration','GroupingVariables',{'Subject','Condition'})
  327. varfun(@numel, tbl,'InputVariables','peak','GroupingVariables',{'Subject','Condition'})
  328. varfun(@numel, tbl,'InputVariables','timeToPeak','GroupingVariables',{'Subject','Condition'})
  330. % % Baseline mixed-effects model - Is there a fixed effect of Condition
  331. % accounting for subject-to-subject differences?
  332. % For each dependent variable, start with a random intercept model
  333. lme1 = fitlme(tbl, 'duration ~ Condition + (1|Subject)')
  335. % If each subject may respond differently to condition manipution,
  336. % include random slopes
  337. lme1_2 = fitlme(tbl, 'duration ~ Condition + (Condition|Subject)')
  338. % Compare models
  339. compare(lme1,lme1_2)
  340. % The small pvalue and large likelihood ratio stat show that including
  341. % random slopes is better, but stick with
  342. anova(lme1_2)
  343. anova(lme1)
  345. % Effect sizes
  346. beta = lme1.Coefficients.Estimate(2)
  347. sd_resid = sqrt(lme1.MSE)
  348. d = beta / sd_resid
  349. beta = lme1_2.Coefficients.Estimate(2)
  350. sd_resid = sqrt(lme1_2.MSE)
  351. d = beta / sd_resid
  353. % % Handling 3 dependent variables
  354. % Separate models + correction
  355. % Fit one model per DV:
  356. lme1 = fitlme(tbl, 'duration ~ Condition + (1|Subject)');
  357. lme2 = fitlme(tbl, 'peak ~ Condition + (1|Subject)');
  358. lme3 = fitlme(tbl, 'timeToPeak ~ Condition + (1|Subject)');
  359. % Correct for multiple testing
  360. p = [lme1.Coefficients.pValue(2), lme2.Coefficients.pValue(2), lme3.Coefficients.pValue(2)];
  361. p_fdr = mafdr(p,'BHFDR',true);
  363. % Visualization
  364. tiledlayout(1,3)
  365. nexttile
  366. boxchart(tbl,"Condition","duration",'MarkerStyle','none','BoxFaceColor','k')
  367. ax=gca; ax.TickDir='out'; ax.LineWidth=1; ax.FontSize=14; box off
  368. ylabel('Duration (s)');
  369. ax.XTickLabel = {'Passive','Intentional'};
  370. xlabel('');
  371. hold on
  372. for s = categories(tbl.Subject)'
  373. idx = tbl.Subject == s{1};
  374. means = groupsummary(tbl(idx,:), 'Condition', 'mean', 'duration');
  375. plot([1 2], means.mean_duration, '-o', 'Color',[.7 .7 .7]);
  376. end
  377. nexttile
  378. boxchart(tbl,"Condition","peak",'MarkerStyle','none','BoxFaceColor','k')
  379. ax=gca; ax.TickDir='out'; ax.LineWidth=1; ax.FontSize=14; box off
  380. ylabel('Peak airflow (au)');
  381. ax.XTickLabel = {'Passive','Intentional'};
  382. xlabel('');
  383. hold on
  384. for s = categories(tbl.Subject)'
  385. idx = tbl.Subject == s{1};
  386. means = groupsummary(tbl(idx,:), 'Condition', 'mean', 'peak');
  387. plot([1 2], means.mean_peak, '-o', 'Color',[.7 .7 .7]);
  388. end
  389. nexttile
  390. boxchart(tbl,"Condition","timeToPeak",'MarkerStyle','none','BoxFaceColor','k')
  391. ax=gca; ax.TickDir='out'; ax.LineWidth=1; ax.FontSize=14; box off
  392. ylabel('Time to peak (ms)');
  393. ax.XTickLabel = {'Passive','Intentional'};
  394. xlabel('');
  395. hold on
  396. for s = categories(tbl.Subject)'
  397. idx = tbl.Subject == s{1};
  398. means = groupsummary(tbl(idx,:), 'Condition', 'mean', 'timeToPeak');
  399. plot([1 2], means.mean_timeToPeak, '-o', 'Color',[.7 .7 .7]);
  400. end
  401. set(gcf,'pos',[300 300 1658 508]);
  403. %% Fig S4C - spectrogram of subset of passive inhales matched for duration of intentional sniffs
  404. % 1b. Bin intentional sniffs and passive inhales by overlapping duration within each subject
  405. % Apart from Subject 2, all subjects show increase in duration for
  406. % intentional vs. passive sniff.
  407. % So try sorting sniff duration for each condition, and take shortest
  408. % duration bin from each subject, except for subject 2 take longest
  409. % duration bin
  411. % Try for 20 subsampled trials of passive inhales for each subject
  412. % Ascending order by default
  413. [inhaleDur_passive_bySub_sorted, inhaleDur_passive_bySub_sortInd] = cellfun(@sort,inhaleDur_passive,'UniformOutput',false);
  415. % Bin sniff onsets by the 20 trials of shortest duration for 5 of 6 subs
  416. sub_inds = [1,3,4,5,6];
  417. for i = 1:length(sub_inds)
  418. sniffonsets_subset_passive(:,sub_inds(i)) = sniffonsets_Passive{sub_inds(i)}(inhaleDur_passive_bySub_sortInd{sub_inds(i)}(1:20));
  419. end
  420. sniffonsets_subset_passive(:,2) = sniffonsets_Passive{2}(inhaleDur_passive_bySub_sortInd{2}(end-19:end));
  422. % Caclulate spectrogram by subjects for the subset of passive inhales
  423. for sub = 1:num_subs
  424. ts = lfp_Passive{sub};
  425. sniffonset_conds = { sniffonsets_subset_passive(:,sub) };
  426. para = [];
  427. para.freq = logspace( log10( 1), log10( 200), 100); % center frequencies for filtering
  428. para.bandwidth = logspace( log10( 2), log10( 30), 100); % bandwidth of each center frequency
  429. para.srate = srate; % sampling rate, in Hz
  430. para.events = sniffonset_conds; % in samples
  431. para.epoch = [-1, 5]; % in seconds, relative to events
  432. para.baseline = [-0.6, -0.1]; % in seconds, relative to events
  433. para.smooth_win = 0.01; % temporal smoothing, in secondsj
  434. para.nbpermuts = 10000; % number of permutations,
  435. para.measure = 'amplitude'; % or 'power'
  436. para.baseline_condition = []; % which condition to use as baseline (empty for just before each of multiple conditions)
  437. para.trl = 'mean'; % 'mean' o r'median' value across trials
  438. [stats_passive{sub}, epoch_time_passive{sub}, ev_epochs_passive{sub}, ~, sd_passive{sub}] = PermutAmpChange_Commonbaseline(ts, para);
  439. end
  440. save('spectro_passive_subsetMatchDur_bySub','para','stats_passive','ev_epochs_passive','epoch_time_passive','sd_passive','-v7.3');
  442. % Combined subs spectrogram for subset of Passive Inhales
  443. % % Subset is for 20 trials that are shortest duration inhale which overlap
  444. % with duration of intentional sniffs, except for sub2 they are the fastest
  446. % Load in previously calculated cued sniff spectrograms to compare
  447. % Should have variable stats_6normosmic_cuedSniff
  449. for sub = 1:6 % 6 normosmics
  450. z_cuedSniff(sub,:,:) = data_supp.stats_6normosmic_cuedSniff{sub}.Z{1};
  451. z_passiveSubset(sub,:,:) = stats_passive{sub}.Z{1};
  452. end
  454. % Average spectrogram output across subjects
  455. z_mean_cuedSniff = squeeze(mean(z_cuedSniff,1));
  456. z_mean_passiveSubset = squeeze(mean(z_passiveSubset,1));
  458. P.P_passive = G_z2p(z_mean_passiveSubset);
  459. FDR.FDR_passive = fdr(P.P_passive(:));
  460. thresh.thresh_passive = G_p2z(FDR.FDR_passive);
  462. P = []; FDR = []; thresh=[];
  463. P.P_cuedSniff = G_z2p(z_mean_cuedSniff);
  464. FDR.FDR_cuedSniff = fdr(P.P_cuedSniff(:));
  465. thresh.thresh_cuedSniff = G_p2z(FDR.FDR_cuedSniff);
  467. for sub = 1:num_subs
  468. resp_epochs_passiveSubset_bySub{sub} = squeeze(G_EpochTS(resp_Passive{sub},sniffonsets_subset_passive(:,sub),[-1,5]*srate));
  469. resp_epochs_cuedSniff_bySub{sub} = squeeze(G_EpochTS(resp_CuedSniff{sub},sniffonsets_CuedSniff{sub},[-1,5]*srate));
  470. end
  471. resp_epochs_passiveSubset_acrossSubs = cell2mat(resp_epochs_passiveSubset_bySub);
  472. resp_mean_epoch_passiveSubset = mean(resp_epochs_passiveSubset_acrossSubs,2);
  473. resp_epochs_cuedSniff_acrossSubs = cell2mat(resp_epochs_cuedSniff_bySub);
  474. resp_mean_epoch_cuedSniff = mean(resp_epochs_cuedSniff_acrossSubs,2);
  476. % Calculate amplitude difference between Passive inhale vs. Cued sniff combined subjects
  477. % Combine ev_epochs across subjects
  478. ev_epochs_6normosmic_cuedSniff = data_supp.ev_epochs_6normosmic_cuedSniff;
  479. ev_epochs_6normosmics_cuedSniff = cat( 3,ev_epochs_6normosmic_cuedSniff{1}{1},ev_epochs_6normosmic_cuedSniff{2}{1},ev_epochs_6normosmic_cuedSniff{3}{1},...
  480. ev_epochs_6normosmic_cuedSniff{4}{1},ev_epochs_6normosmic_cuedSniff{5}{1},ev_epochs_6normosmic_cuedSniff{6}{1} );
  481. ev_epochs_6normosmics_passiveSubset = cat( 3,ev_epochs_passive{1}{1},ev_epochs_passive{2}{1},ev_epochs_passive{3}{1},...
  482. ev_epochs_passive{4}{1},ev_epochs_passive{5}{1},ev_epochs_passive{6}{1} );
  484. tic
  485. mat = cat(3,ev_epochs_6normosmics_cuedSniff,ev_epochs_6normosmics_passiveSubset);
  486. mat = permute(mat, [3 1 2]);
  487. g = [ones(size(ev_epochs_6normosmics_cuedSniff,3),1) ; 2*ones(size(ev_epochs_6normosmics_passiveSubset,3),1)];
  488. AmpChange_cuedVsPassiveSubsetMatchedDur = G_PowChange_Clust(mat, g, epoch_time_passive{1}, para.freq, 'baseline', para.baseline,'db_transf', 'no','common_bs','yes');
  489. save('AmpChange_cuedVsPassiveSubsetMatchedDur','AmpChange_cuedVsPassiveSubsetMatchedDur','para','-v7.3');
  490. toc
  492. AmpChange_P = []; AmpChange_FDR=[]; AmpChange_thresh=[];
  493. AmpChange_P = G_z2p(AmpChange_cuedVsPassiveSubsetMatchedDur.statsmap);
  494. AmpChange_FDR = fdr(AmpChange_P(:));
  495. AmpChange_thresh = G_p2z( AmpChange_FDR)
  497. % Plot mean spectrogram for Cued Sniff vs. Passive Inhale averaged across subjects
  498. Time = -1 : 1/srate : 5;
  499. tiledlayout(1,3)
  500. nexttile
  501. [h, c, clusts, clustinfo]=G_Imagesc(z_mean_cuedSniff,'yaxis',freq,...
  502. 'xaxis', Time, 'xtick', [-0.3,0:0.5:5],'threshold',thresh.thresh_cuedSniff,...
  503. 'ytick', [1,2,3,8,15,30,50,100,200],'Fontsize',16,'overlay',resp_mean_epoch_cuedSniff,...
  504. 'clim',[-3,3],'ctick', [-3,0,3],...
  505. 'ytick_precision',0,'xlabel','Time (s)','ylabel','Frequency (Hz)');
  506. ax=gca;
  507. ax.FontSize = 16; ax.LineWidth = 1; box off
  508. xlim([0.69*srate,3.51*srate]);
  509. title('Intentional Sniff');
  510. colorbar off
  511. nexttile
  512. G_Imagesc(z_mean_passiveSubset,'yaxis',freq,...
  513. 'xaxis', Time, 'xtick', [-0.3,0:0.5:5],'threshold',thresh.thresh_passive,...
  514. 'ytick', [1,2,3,8,15,30,50,100,200],'Fontsize',16,'overlay',resp_mean_epoch_passiveSubset,'overlay_ref',resp_mean_epoch_cuedSniff,...
  515. 'clim',[-3,3],'ctick', [-3,0,3],'clabel','Amplitude (z score)',...
  516. 'ytick_precision',0);
  517. xlim([0.69*srate,3.51*srate]);
  518. title('Passive Inhale');
  519. ylabel('');
  520. xlabel('');
  521. ax=gca;
  522. ax.FontSize = 16; ax.LineWidth = 1; box off
  523. % set(gcf,'pos',[758 784 1047 478]); % without difference plot
  524. nexttile
  525. h=G_Imagesc(AmpChange_cuedVsPassiveSubsetMatchedDur.statsmap,'yaxis',para.freq,...
  526. 'xaxis', Time, 'xtick', [-0.3,0:0.5:5], ...
  527. 'ytick', [1,2,8,15,30,50,100,200],'Fontsize',16,...
  528. 'clim',[-3,3],'ctick', [-3,0,3],'clabel','Amplitude Difference (z score)',...
  529. 'ytick_precision',0);
  530. xlabel(''); ylabel('');
  531. title('Difference');
  532. for k = 1:length(AmpChange_cuedVsPassiveSubsetMatchedDur.clust{1}.both.info.PixelIdxList)
  533. mask = zeros(AmpChange_cuedVsPassiveSubsetMatchedDur.clust{1}.both.info.ImageSize);
  534. if AmpChange_cuedVsPassiveSubsetMatchedDur.clust{1}.both.statssum_pval(k) < 0.05
  535. mask(AmpChange_cuedVsPassiveSubsetMatchedDur.clust{1}.both.info.PixelIdxList{k}) = 1;
  536. G_Outline_Pixels(h,mask,'linecolor','k','LineWidth',1.5);
  537. end
  538. end
  539. xlim([0.69*srate,3.51*srate]);
  540. set(gcf,'pos',[758 784 1499 478]); % With difference plot
  541. ax.FontSize = 16; ax.LineWidth = 1; box off
  543. %% Fig S4D to E
  545. ph_epochs_byRespPhaseBin = data_supp.ph_epochs_byRespPhaseBin;
  546. ph_epochs_byRespPhase_combinedSubs{1} = cat(3,ph_epochs_byRespPhaseBin{1}{1}(:,:,inhaleDur_passive_bySub_sortInd{1}),ph_epochs_byRespPhaseBin{1}{2}(:,:,inhaleDur_passive_bySub_sortInd{2}),...
  547. ph_epochs_byRespPhaseBin{1}{3}(:,:,inhaleDur_passive_bySub_sortInd{3}),ph_epochs_byRespPhaseBin{1}{4}(:,:,inhaleDur_passive_bySub_sortInd{4}),...
  548. ph_epochs_byRespPhaseBin{1}{5}(:,:,inhaleDur_passive_bySub_sortInd{5}),ph_epochs_byRespPhaseBin{1}{6}(:,:,inhaleDur_passive_bySub_sortInd{6}));
  550. ph_epochs_byRespPhase_combinedSubs{2} = cat(3,ph_epochs_byRespPhaseBin{2}{1},ph_epochs_byRespPhaseBin{2}{2},ph_epochs_byRespPhaseBin{2}{3},...
  551. ph_epochs_byRespPhaseBin{2}{4},ph_epochs_byRespPhaseBin{2}{5},ph_epochs_byRespPhaseBin{2}{6});
  553. ph_epochs_cuedSniff = data_supp.ph_epochs_byRespPhase_combinedSubs{2};
  554. ph_epochs_passive = cat(3,ph_epochs_byRespPhaseBin{1}{1},ph_epochs_byRespPhaseBin{1}{2},...
  555. ph_epochs_byRespPhaseBin{1}{3},ph_epochs_byRespPhaseBin{1}{4},...
  556. ph_epochs_byRespPhaseBin{1}{5},ph_epochs_byRespPhaseBin{1}{6});
  558. param.freq = logspace( log10( 1), log10( 100), 80); % logarithmic steps in frequency for filtering
  559. param.bandwidth = logspace( log10( 2), log10( 20), 80); % logarithmic steps in bandwidth for filtering
  561. % Calculate ITPC difference between combined subs respiratory-normalized conditions
  562. % % Calculate ITPC difference between Cued sniff and Passive
  563. ITPCdiff_normResp_combinedSubs = data_supp.ITPC_normResp_combinedSubs{2} - data_supp.ITPC_normResp_combinedSubs{1};
  564. ITPCdiff_fisherZ_normResp_combinedSubs = atanh(ITPCdiff_normResp_combinedSubs);
  566. % % Permutation of ITPC over circularly shifted data, 1000 iterations
  567. numpts = size(ph_epochs_byRespPhase_combinedSubs{2},2);
  568. num_perms = 1000;
  569. for perm = 1:num_perms
  570. fprintf( 'Permutation: %d/%d \r', perm, num_perms);
  571. ph_tmp_passive=ph_epochs_byRespPhase_combinedSubs{1};
  572. ph_tmp_cuedSniff=ph_epochs_byRespPhase_combinedSubs{2};
  573. for trial = 1:size(ph_epochs_byRespPhase_combinedSubs{1},3)
  574. shift_amount = circshift( 1:numpts, randsample( 1:numpts, 1));
  575. ph_tmp_passive( :, :, trial) = ph_epochs_byRespPhase_combinedSubs{1}( :, shift_amount, trial );
  576. end
  577. for trial = 1:size(data_supp.ph_epochs_byRespPhase_combinedSubs{2},3)
  578. shift_amount = circshift( 1:numpts, randsample( 1:numpts, 1));
  579. ph_tmp_cuedSniff( :, :, trial) = ph_epochs_byRespPhase_combinedSubs{2}( :, shift_amount, trial);
  580. end
  581. permuteITPC_passive_tmp = []; permuteITPC_cuedSniff_tmp = [];
  582. permuteITPC_passive_tmp(:,:,perm) = squeeze(abs(mean(exp(1i * ph_tmp_passive), 3)));
  583. permuteITPC_cuedSniff_tmp(:,:,perm) = squeeze(abs(mean(exp(1i * ph_tmp_cuedSniff), 3)));
  584. permuteITPCdiff_normResp_combinedSubs(:,:,perm) = permuteITPC_cuedSniff_tmp(:,:,perm) - permuteITPC_passive_tmp(:,:,perm);
  585. permuteITPCdiff_fisherZ_passiveVcuedsniff_normResp(:,:,perm) = atanh(permuteITPCdiff_normResp_combinedSubs(:,:,perm));
  586. end
  588. % % % Cluster-based correction
  589. [h_clust_normResp_combinedSubs_ITPCdiff, p_clust_normResp_combinedSubs_ITPCdiff, clustinfo_normResp_combinedSubs_ITPCdiff] = cluster_test(ITPCdiff_fisherZ_normResp_combinedSubs, permuteITPCdiff_fisherZ_passiveVcuedsniff_normResp);
  591. % Plot ITPC difference between Cued Sniff and passive with Resp Overlay
  593. % Respiratory epochs for overlay
  594. resp_epochs_byRespPhaseBin = data_supp.resp_epochs_byRespPhaseBin;
  595. resp_epochs_byRespPhaseBin_allsubs_cuedSniff = squeeze(cat(3,resp_epochs_byRespPhaseBin{2}{1},resp_epochs_byRespPhaseBin{2}{2},resp_epochs_byRespPhaseBin{2}{3},...
  596. resp_epochs_byRespPhaseBin{2}{4},resp_epochs_byRespPhaseBin{2}{5},resp_epochs_byRespPhaseBin{2}{6}));
  597. resp_meanEpoch_byRespPhaseBin_cuedSniff = mean(resp_epochs_byRespPhaseBin_allsubs_cuedSniff,2);
  599. resp_epochs_byRespPhaseBin_allsubs_Passive = squeeze(cat(3,resp_epochs_byRespPhaseBin{1}{1},resp_epochs_byRespPhaseBin{1}{2},resp_epochs_byRespPhaseBin{1}{3},...
  600. resp_epochs_byRespPhaseBin{1}{4},resp_epochs_byRespPhaseBin{1}{5},resp_epochs_byRespPhaseBin{1}{6}));
  601. resp_meanEpoch_byRespPhaseBin_Passive = mean(resp_epochs_byRespPhaseBin_allsubs_Passive,2);
  603. [redblue,red,blue]=CustomColormap;
  604. % subplot(3,2,sub);
  605. h=G_Imagesc(ITPCdiff_fisherZ_normResp_combinedSubs,'yaxis',param.freq,...
  606. 'ytick',[1,2,3,8,15,30,50,100],'xtick',1:50:202,'ytick_precision',1,'clim',[0,0.5],...
  607. 'ylabel','Frequency (Hz)','xlabel','','overlay',resp_meanEpoch_byRespPhaseBin_cuedSniff);
  608. title(sub_IDs{sub});
  609. for k = 1:length(clustinfo_normResp_combinedSubs_ITPCdiff.pos_clusters)
  610. if clustinfo_normResp_combinedSubs_ITPCdiff.pos_clusters(k).p < 0.05
  611. G_Outline_Pixels(h,clustinfo_normResp_combinedSubs_ITPCdiff.pos_clusters(k).inds,'linecolor','k','linewidth',1.5);
  612. end
  613. end
  614. colormap(red);
  615. title('Intentional Sniff vs. Passive Nasal Inhalation');
  616. ax=gca;
  617. ax.XTickLabel={'-\pi','','0','','\pi'};
  618. ax.FontSize = 16;
  619. ax.XLabel.String='Respiratory Phase';
  620. cb=colorbar;
  621. cb.Label.String='ITPC difference value (Fisher-transformed)';
  622. cb.Limits = [0,0.5];
  623. cb.Ticks = 0:0.25:1;
  624. % cb.Label.String='ITPCz';
  625. % cb.Limits = [0,16];
  627. % %% Plot combined subjects cued sniff and passive normalized by respiratory phase
  628. [redblue,red,blue]=CustomColormap;
  629. tiledlayout(1,2)
  630. nexttile
  631. h=G_Imagesc(ITPC_normResp_combinedSubs{2},'yaxis',param.freq,...
  632. 'ytick',[1,2,3,8,15,30,50,100],'xtick',1:50:202,'ytick_precision',1,'clim',[0,0.5],...
  633. 'ylabel','Frequency (Hz)','xlabel','','overlay',resp_meanEpoch_byRespPhaseBin_cuedSniff);
  634. title(sub_IDs{sub});
  635. for k = 1:length(clustinfo_normResp_combinedSubs{2}.pos_clusters)
  636. if clustinfo_normResp_combinedSubs{2}.pos_clusters(k).p < 0.05
  637. G_Outline_Pixels(h,clustinfo_normResp_combinedSubs{2}.pos_clusters(k).inds,'linecolor','k','linewidth',1.5);
  638. end
  639. end
  640. colormap(red);
  641. title('Intentional Sniff');
  642. ax=gca;
  643. ax.XTickLabel={'-\pi','','0','','\pi'};
  644. ax.FontSize = 16;
  645. ax.XLabel.String='Respiratory Phase';
  646. cb=colorbar;
  647. cb.Label.String='ITPC';
  648. cb.Limits = [0,0.5];
  649. cb.Ticks = 0:0.25:1;
  650. colorbar off
  651. % Plot combined subjects passive normalized by respiratory phase
  652. [redblue,red,blue]=CustomColormap;
  653. nexttile
  654. % subplot(3,2,sub);
  655. h=G_Imagesc(ITPC_normResp_combinedSubs{1},'yaxis',param.freq,...
  656. 'ytick',[1,2,3,8,15,30,50,100],'xtick',1:50:202,'ytick_precision',1,'clim',[0,0.5],...
  657. 'ylabel','Frequency (Hz)','xlabel','','overlay',resp_meanEpoch_byRespPhaseBin_Passive);
  658. title(sub_IDs{sub});
  659. for k = 1:length(clustinfo_normResp_combinedSubs{1}.pos_clusters)
  660. if clustinfo_normResp_combinedSubs{1}.pos_clusters(k).p < 0.05
  661. G_Outline_Pixels(h,clustinfo_normResp_combinedSubs{1}.pos_clusters(k).inds,'linecolor','k','linewidth',1.5);
  662. end
  663. end
  664. colormap(red);
  665. title('Passive Inhalation');
  666. ax=gca;
  667. ax.XTickLabel={'-\pi','','0','','\pi'};
  668. ax.FontSize = 16;
  669. ax.XLabel.String='Respiratory Phase';
  670. cb=colorbar;
  671. cb.Label.String='ITPC';
  672. cb.Limits = [0,0.5];
  673. cb.Ticks = 0:0.25:1;
  674. xlabel('');
  675. ylabel('');
  676. % cb.Label.String='ITPCz';
  677. % cb.Limits = [0,16];
  678. % saveas(gcf,'ITPC_Passive_normResp_combined6subs','fig');
  679. % saveas(gcf,'ITPC_Passive_normResp_combined6subs','png');
  680. set(gcf,'pos',[537 818 1023 420]);
  682. %% ITPC Permuted Distribution for Cued sniff vs. Passive inhale Subset matched inhale dur
  683. total_trials_cued = size(ph_epochs_cuedSniff,3);
  684. total_trials_passive = size(ph_epochs_passive,3);
  685. freq = logspace( log10( 1), log10( 100), 80); % logarithmic steps in frequency for filtering
  686. theta_ind = find(freq>=2&freq<=8);
  687. tic
  688. for i = 1:1000
  689. % rand_passive_inds = randperm( total_trials, num_fixedTrials ); % old
  690. rand_inds_cued = randi( total_trials_cued, [1, total_trials_cued] );
  691. ph_epochs_resample_cuedSniff = ph_epochs_cuedSniff(:,:,rand_inds_cued);
  692. ITPC_resample_cuedSniff{i} = squeeze(abs(mean(exp(1i * ph_epochs_resample_cuedSniff), 3))); % length of mean vector over trials (dimension 3 of phase mat) = ITPC
  693. ITPC_distr_maxTheta_cuedSniff(i) = max(ITPC_resample_cuedSniff{i}(theta_ind,:),[],'all');
  694. rand_inds_passive = randi( total_trials_passive, [1, total_trials_passive] );
  695. ph_epochs_resample_passive = ph_epochs_passive(:,:,rand_inds_passive);
  696. ITPC_resample_passive{i} = squeeze(abs(mean(exp(1i * ph_epochs_resample_passive), 3))); % length of mean vector over trials (dimension 3 of phase mat) = ITPC
  697. ITPC_distr_maxTheta_passive(i) = max(ITPC_resample_passive{i}(theta_ind,:),[],'all');
  698. end
  699. toc
  701. % Plot resampled distribution of passive vs. intentional sniff
  702. histogram(ITPC_distr_maxTheta_cuedSniff,32,'FaceColor',[0 0.6 1])
  703. hold on
  704. histogram(ITPC_distr_maxTheta_passive,32,'FaceColor',[1 0.6 0])
  705. legend('Intentional Sniff', 'Passive Inhale','box','off');
  706. ax=gca; ax.TickDir = 'out'; ax.LineWidth=1; box off; ax.FontSize=16;
  707. xlabel('ITPC');
  708. ylabel('Resample #');
  709. saveas(gcf,'ITPC_resampled_cuedVsPassiveSubsetMatchedDur','fig');
  710. exportgraphics(gcf,'ITPC_resampled_cuedVsPassiveSubsetMatchedDur.pdf','ContentType','vector');
  712. ITPC_distr_maxTheta_cuedSniff - ITPC_distr_maxTheta_passive;
  713. histogram(ans,32)
  715. % Stats on bootstrapped ITPC difference between Intentional sniff and
  716. % Passive Inhale
  717. bootstrap_ITPCdiff = ITPC_distr_maxTheta_cuedSniff - ITPC_distr_maxTheta_passive;
  719. [muHat_ITPCdiff,sigmaHat_ITPCdiff,muCI_ITPCdiff,sigmaCI_ITPCdiff] = normfit(bootstrap_ITPCdiff, 0.05)
  721. %% Fig S4F
  722. ev_epochs_6normosmic_Passive = data_supp.ev_epochs_6normosmic_Passive;
  723. sd_6_normosmic_Passive = data_supp.sd_6_normosmic_Passive;
  724. para = data_supp.para;
  726. for sub = 1:6 % first 6 subs are normosmic
  727. resp_epochs_Passive_bySub{sub} = squeeze(zscore(G_EpochTS(data_supp.resp_Passive{sub}, (data_supp.sniffonsets_Passive{sub}, [-1,5]*srate)));
  728. resp_meanEpoch_Passive_bySub{sub} = squeeze(mean(resp_epochs_Passive_bySub{sub},2));
  729. end
  731. % Passive - first 6 subs are normosmics
  732. baseline = para.baseline;
  733. baselineloc = G_RangeLoc( epoch_time{1}, baseline);
  734. freq = para.freq;
  735. % theta_ind = find(freq >= 3 & freq <= 8);
  736. theta_ind = find(freq > 2 & freq <= 8);
  737. num_freqs = size(ev_epochs_6normosmic_Passive{1}{1},1);
  738. numtrials=[]; avg_Passive=[]; response_Passive=[]; amp_epochs_corec_Passive=[]; amp_epochs_theta_Passive=[]; mean_amp_theta_Passive=[]; sem_amp_theta_Passive=[];
  739. for sub = 1:6 % first 6 subs are normosmics
  740. numtrials(sub) = size(ev_epochs_6normosmic_Passive{sub}{1},3);
  741. avg_Passive{sub} = mean(ev_epochs_6normosmic_Passive{sub}{1}(:,baselineloc),2);
  742. rand_passive_inds_sub{sub} = randperm( size(ev_epochs_6normosmic_Passive{sub}{1},3),size(ev_epochs_6normosmic_cuedSniff{sub}{1},3) );
  743. response_Passive{sub} = bsxfun( @minus, ev_epochs_6normosmic_Passive{sub}{1}, avg_Passive{sub} );
  744. response_theta_Passive{sub} = squeeze(mean(response_Passive{sub}(theta_ind,:,:),1));
  745. amp_epochs_corec_Passive{sub} = bsxfun( @rdivide, response_Passive{sub}, sd_6_normosmic_Passive{sub}{1} );
  746. amp_epochs_theta_Passive{sub} = squeeze(mean(amp_epochs_corec_Passive{sub}(theta_ind,:,:),1));
  747. inc_theta_Passive{sub} = squeeze(mean(amp_epochs_theta_Passive{sub}(post500ms_ind,rand_passive_inds_sub{sub}) - amp_epochs_theta_Passive{sub}(pre500ms_ind,rand_passive_inds_sub{sub}),1));
  748. inc_theta_Passive{sub} = mean(inc_theta_Passive{sub},1);
  749. mean_amp_theta_Passive(:,sub) = mean(amp_epochs_theta_Passive{sub},2);
  750. sem_amp_theta_Passive(:,sub) = std(amp_epochs_theta_Passive{sub},[],2) / sqrt(numtrials(sub));
  751. end
  753. % collapse across subjects
  754. thetaAmp_epochs_Passive=[]; resp_epochs_Passive_acrossSubs=[];
  755. thetaAmp_epochs_Passive = cell2mat(amp_epochs_theta_Passive);
  756. resp_epochs_Passive_acrossSubs = cell2mat(resp_epochs_Passive_bySub(1:6));
  757. inc_theta_Passive_allSubs = cell2mat(inc_theta_Passive);
  758. inc_theta_Passive_allSubs = cell2mat(inc_theta_Passive);
  760. % calculate inhale peak, latency to inhale peak, theta peak,
  761. % latency to theta peak for Passive
  762. inhalePk_Passive=[]; time2inhalePk_ind_Passive=[]; thetaPk_CuedSniff=[]; time2thetaPk_ind_Passive=[]; numTrials_Passive_allSubs=[];
  763. numTrials_Passive_allSubs = size(thetaAmp_epochs_Passive,2);
  764. for trial = 1:size(resp_epochs_Passive_acrossSubs,2)
  765. [inhalePk_Passive(trial), time2inhalePk_ind_Passive(trial)] = max(resp_epochs_Passive_acrossSubs(500:1500,trial));
  766. [thetaPk_Passive(trial), time2thetaPk_ind_Passive(trial)] = max(thetaAmp_epochs_Passive(500:1500,trial));
  767. end
  768. time2inhalePk_ms_Passive=[]; time2thetaPk_ms_Passive=[];
  769. time2inhalePk_ms_Passive = ((time2inhalePk_ind_Passive) / srate) * 1000;
  770. time2thetaPk_ms_Passive = ((time2thetaPk_ind_Passive) / srate) * 1000;
  772. isoutlier(thetaPk_Passive);
  773. thetaPk_Passive_outliersRem = thetaPk_Passive(~isoutlier(thetaPk_Passive));
  774. inhalePk_Passive_outliersRem = inhalePk_Passive(~isoutlier(thetaPk_Passive));
  776. [r,p] = corr(thetaPk_Passive_outliersRem',inhalePk_Passive_outliersRem')
  778. scatter(inhalePk_Passive_outliersRem', thetaPk_Passive_outliersRem', 'filled','MarkerEdgeColor', [0 0 0.7],'MarkerFaceColor',[0 0 0.7])
  779. ax=gca; ax.TickDir='out'; ax.LineWidth = 1; box off; ax.FontSize=16;
  780. xlabel('Peak inhale magnitude (z score)');
  781. ylabel('Maximal Theta Amplitude (normalized)');
  782. ax = gca; ax.TickDir = 'out'; ax.FontSize = 16; ax.LineWidth = 1;
  783. text(0.75,27.5,append('r = ',num2str(r,2)),'FontSize',14)
  784. text(0.75,25,append('p = ',num2str(p,1)),'FontSize',14)
  785. l=lsline;
  786. l.LineWidth = 2;
  787. set(gcf,'pos',[1000 728 815 510]);
  789. %% Fig S5 - Modulation index on proximal vs. distal contacts
  790. % Filtered steps for phase frequency across theta band - ODOR
  791. para.ph_freq = linspace(2,8,19);
  792. para.ph_bandwidth = 2;
  793. para.amp_freq = logspace(log10(29),log10(153),32);
  794. para.amp_bandwidth = logspace(log10(15),log10(42),32);
  795. nbamps = length( para.amp_freq);
  796. nbphs = length(para.ph_freq);
  797. ph=[]; amp=[]; ph_epochs=[]; amp_epochs=[]; ph_concat=[]; amp_concat=[];
  798. ph_allTrials=[]; amp_allTrials=[]; amp_ts=[]; ph_ts=[];
  799. for sub = 1:num_subs
  800. for freq_idx = 1 : nbphs
  801. bpfreq = para.ph_freq( freq_idx) + [-0.5, 0.5]*para.ph_bandwidth;
  802. % fprintf( 'phase filtering %d %d\n', bpfreq(1), bpfreq(2));
  803. if bpfreq(1) < eps
  804. ph_ts{sub}( freq_idx, :) = ft_preproc_lowpassfilter( distalCh_odorTasks{sub}, srate, bpfreq(2), [], 'fir');
  805. else
  806. ph_ts{sub}( freq_idx, :) = ft_preproc_bandpassfilter( distalCh_odorTasks{sub}, srate, bpfreq, [], 'fir');
  807. end
  808. end
  809. ph{sub} = angle( hilbert( ph_ts{sub}'));
  810. ph{sub} = wrapTo2Pi( ph{sub});
  811. ph{sub} = ph{sub}';
  812. ph_epochs{sub} = squeeze(G_EpochTS(ph{sub}, sniffonsets_odorTasks{sub}, [1,4]*srate));
  813. % ph_epochs{sub} = squeeze(G_EpochTS(ph{sub}, trough_thresh_times{sub}, [-0.25,0.25]*srate));
  815. for freq_idx = 1 : nbamps
  816. bpfreq = para.amp_freq( freq_idx) + [-0.5, 0.5]*para.amp_bandwidth( freq_idx);
  817. fprintf( 'amplitude filtering %d %d\n', bpfreq(1), bpfreq(2));
  818. amp_ts{sub}( freq_idx, :) = ft_preproc_bandpassfilter( lfp_odorTasks{sub}, srate, bpfreq, [], 'fir');
  819. end
  820. amp{sub} = abs(hilbert(amp_ts{sub}));
  821. amp_epochs{sub} = squeeze(G_EpochTS(amp{sub}, sniffonsets_odorTasks{sub}, [1,4]*srate));
  822. % amp_epochs{sub} = squeeze(G_EpochTS(amp{sub}, trough_thresh_times{sub}, [-0.25,0.25]*srate));
  823. end
  825. % Concatenate phase and amplitude epochs
  826. ph_allTrials = cat(3,ph_epochs{1},ph_epochs{2},ph_epochs{3},ph_epochs{4},...
  827. ph_epochs{5}, ph_epochs{6});
  828. ph_concat = reshape(ph_allTrials, size(ph_allTrials,1), []);
  829. amp_allTrials = cat(3,amp_epochs{1},amp_epochs{2},amp_epochs{3},...
  830. amp_epochs{4},amp_epochs{5},amp_epochs{6});
  831. amp_concat = reshape(amp_allTrials,length(para.amp_freq),[]);
  833. % Calculate MI
  834. para.nbsurrogate = 200;
  835. para.nbbins = 18;
  836. [raw_MI_distal, norm_MI_distal, bin_amp_distal, ph_wins_distal] = G_MI( ph_concat, amp_concat, para.nbsurrogate, para.nbbins, srate);
  838. % Plot proximal vs. distal
  839. norm_MI_proximal = data_supp.MI_proximal.norm_MI;
  840. norm_MI_p = G_z2p(norm_MI_proximal);
  841. norm_MI_FDR = fdr(norm_MI_p(:));
  842. norm_MI_thresh = G_p2z(norm_MI_FDR)
  844. norm_MI_distal = norm_MI;
  845. plot(para.amp_freq,norm_MI_proximal(:,7),'k','LineWidth',2)
  846. yline(norm_MI_thresh(2),'Color',[0.7 0 0],'LineStyle','--','LineWidth',1)
  847. hold on
  848. plot(para.amp_freq,norm_MI_distal(:,7),'Color',[0.6 0.6 0.6],'LineWidth',2)
  849. legend({'Proximal','','Distal'},'box','off');
  850. ax=gca; ax.FontSize = 16; ax.LineWidth=1; ax.TickDir='out'; box off
  851. xlim([29 150]);
  852. xlabel('Frequency (Hz)');
  853. ylabel('Modulation Index (z score)');
  854. set(gcf,'pos',[300 300 673 571]);
  856. %% Fig S6 - Sniffs with and without intranasal electrode in place
  857. %% Fig S6A - Sniff traces
  858. Time = -0.5 : 1/srate : 3;
  860. figure
  861. for sub = 1:7
  862. plot(Time,data_supp.resp_epochs{sub},'Color',[1 0 0 0.1],'LineWidth',1)
  863. hold on
  864. end
  865. for sub = 8:13
  866. plot(Time,data_supp.resp_epochs{sub-7},'Color',[0 0 1 0.1],'LineWidth',1);
  867. hold on
  868. end
  869. % Create "dummy" lines with max alpha (alpha = 1) for the legend only
  870. % Set 'HandleVisibility' to 'off' so they don't show up in the main plot
  871. h1_legend = plot(NaN, NaN, 'r-', 'LineWidth', 2, 'HandleVisibility', 'off');
  872. h2_legend = plot(NaN, NaN, 'b-', 'LineWidth', 2, 'HandleVisibility', 'off');
  873. % Note: The 'HandleVisibility' does not affect how they appear in the legend
  874. % once you explicitly include their handles.
  875. % Create the legend using the handles of the dummy lines
  876. legend([h1_legend, h2_legend], {'Without', 'With'}, 'Location', 'best','box','off');
  877. hold off;
  878. set(gcf,'pos',[300 300 913 732]);
  879. ax=gca; ax.TickDir = 'out'; ax.LineWidth=1; box off; ax.FontSize=16;
  880. ylabel('Airflow (z score)');
  881. xlabel('Time (s)');
  883. %% Fig S6B - histogram of sniff duration with and without intranasal electrode
  884. % Mark exhale onsets to find duration
  886. nbsamps = size(data_supp.resp_epochs{1},1);
  887. base_idx = zeros([1,nbsamps]);
  888. % Time = -1 : 1/srate : 5;
  889. % base_idx(Time >= -1 & Time <= 0) = ones([1,length(0:1*srate)]);
  890. Time = -0.5 : 1/srate : 3.5;
  891. base_idx(1:0.5*srate+1) = ones([1,length(1:0.5*srate+1)]);
  892. resp_meanCorec_epochs = [];
  893. % baseline correct with scale - values fall from 0 to 1
  894. for sub = 1:num_total_subs
  895. resp_bc_epochs{sub} = G_BaselineCorrect(data_supp.resp_epochs{sub},1,base_idx,'scale');
  896. resp_meanCorec_epochs(:,sub) = mean(resp_bc_epochs{sub},2);
  897. end
  900. % Inputs:
  901. % resp_data : 1x13 cell array
  902. % each cell = [samples x trials] matrix
  903. % time : [samples x 1] time vector from -0.5 to 3 s
  905. nSub = numel(resp_bc_epochs);
  907. % Preallocate output
  908. exhale_sample = cell(nSub,1);
  909. exhale_time = cell(nSub,1);
  911. % Define detection window
  912. win_idx = find(Time >= 0.25 & Time <= 3);
  914. for s = 1:nSub
  916. data = resp_bc_epochs{s}; % samples x trials
  917. [nSamp, nTrials] = size(data);
  919. exhale_sample{s} = nan(1,nTrials);
  920. exhale_time{s} = nan(1,nTrials);
  922. for tr = 1:nTrials
  924. y = data(:,tr);
  926. % Restrict to detection window
  927. y_win = y(win_idx);
  928. t_win = Time(win_idx);
  930. % Find falling crossings of 0.5
  931. cross_idx = find(y_win(1:end-1) >= 0.5 & ...
  932. y_win(2:end) < 0.5);
  934. if ~isempty(cross_idx)
  936. % Take first falling crossing
  937. idx = cross_idx(1);
  939. % Linear interpolation for better precision
  940. y1 = y_win(idx);
  941. y2 = y_win(idx+1);
  942. t1 = t_win(idx);
  943. t2 = t_win(idx+1);
  945. frac = (0.5 - y1) / (y2 - y1);
  946. t_cross = t1 + frac * (t2 - t1);
  948. % Store
  949. exhale_time{s}(tr) = t_cross;
  951. % Convert to nearest global sample index
  952. [~, global_idx] = min(abs(Time - t_cross));
  953. exhale_sample{s}(tr) = global_idx;
  954. end
  956. end
  957. end
  959. % Calculate sniff duration for ECoG patients vs. OBE participants
  960. ECoG_IDs = data_supp.resp_epochs_IDs(1:7);
  961. OB_IDs = data_supp.resp_epochs_IDs(8:13);
  962. num_ECoG = length(ECoG_IDs);
  963. num_OB = length(OB_IDs);
  965. for sub = 1:num_ECoG
  966. sniff_duration_ECoG{sub} = (exhale_sample{sub} - 0.5*srate)/srate;
  967. end
  968. meanBySub_sniff_duration_ECoG = cellfun(@mean, sniff_duration_ECoG);
  969. mean_sniff_duration_ECoG = mean(meanBySub_sniff_duration_ECoG)
  970. sem_sniff_duration_ECoG = std(meanBySub_sniff_duration_ECoG,[])/sqrt(num_ECoG);
  972. for sub = 1:num_OB
  973. sniff_duration_OB{sub} = (exhale_sample{sub+num_ECoG} - 0.5*srate) / srate;
  974. end
  975. meanBySub_sniff_duration_OB = cellfun(@mean, sniff_duration_OB);
  976. mean_sniff_duration_OB = mean(meanBySub_sniff_duration_OB);
  977. sem_sniff_duration_OB = std(meanBySub_sniff_duration_OB,[])/sqrt(num_OB);
  979. % Bar plot of sniff duration between OBE and ECoG
  980. [h,p,ci,stats]=ttest2(meanBySub_sniff_duration_OB, meanBySub_sniff_duration_ECoG)
  982. mean_sniffDur_OB_ECoG(1) = mean_sniff_duration_OB;
  983. mean_sniffDur_OB_ECoG(2) = mean_sniff_duration_ECoG;
  984. sem_sniffDur_OB_ECoG(1) = sem_sniff_duration_OB;
  985. sem_sniffDur_OB_ECoG(2) = sem_sniff_duration_ECoG;
  987. % plot_IDs2 = {'','P1','P2','P3','P4','P5','P6'};
  988. bar( mean_sniffDur_OB_ECoG, 'FaceColor','none',...
  989. 'EdgeColor','k','LineWidth',1)
  990. hold on
  991. scatter(1, meanBySub_sniff_duration_OB, 75, sub_colors,'LineWidth',1)
  992. scatter(2, meanBySub_sniff_duration_ECoG, 75, 'k', 'LineWidth',1)
  993. er = errorbar(mean_sniffDur_OB_ECoG, sem_sniffDur_OB_ECoG);
  994. er.Color = [0 0 0];
  995. er.LineStyle = 'none';
  996. er.LineWidth = 1;
  997. set(gcf,'pos',[300 300 566 603]);
  998. xticklabels({'With','Without'});
  999. xlabel('Presence of intranasal electrode');
  1000. ylabel('Sniff duration (s)');
  1001. ax=gca; ax.TickDir='out'; ax.LineWidth=1; ax.FontSize=16; box off

code_suppFigs.m, under CC-BY-4.0 · at the source

Overview

Authors: Andrew Sheriff1, Mahmoud Omidbeigi1, Gregory Lane1, Qiaohan Yang1, Guangyu Zhou1, Adam Dede1, Naelly Arriaga1, Vivek Sagar1, Ania M. Holubecki1, Rodrigo M. Braga1, Leslie M. Kay2,3, Bruce K. Tan4, Christina Zelano1
ORCID iDs: Andrew Sheriff, Gregory Lane, Guangyu Zhou, Adam Dede, Vivek Sagar, Ania M. Holubecki, Rodrigo M. Braga, Leslie M. Kay, Bruce K. Tan, Christina Zelano
  1. Department of Neurology, Northwestern University Feinberg School of Medicine, Chicago, IL, USA
  2. Department of Psychology, The University of Chicago, Chicago, IL, USA
  3. Institute for Mind and Biology, The University of Chicago, Chicago, IL, USA
  4. Department of Otolaryngology—Head and Neck Surgery, Northwestern University Feinberg School of Medicine, Chicago, IL, USA
Institutions: Northwestern University (United States); University of Chicago (United States)
Journal: Science advances, volume 12, issue 27, article eaee1002
Dates: received 20 November 2025; accepted 14 May 2026; published online 3 July 2026; in print July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1126/sciadv.aee1002 · PMID 42397924 · PMCID PMC13330860 · OpenAlex W7167286658
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism)
Methods: Spectral & time-frequency, Preprocessing, Connectivity, Statistics, Single-unit activity, calcium imaging, Smoothing, state filtering, decompositions, Physiology & signal measures
MeSH: Odorants*, Olfactory Bulb*, Olfactory Perception*, Smell*, Theta Rhythm*, Humans (* major topic)
Journal subjects: Neuroscience, Neurophysiology
Topic: Olfactory and Sensory Function Studies (Sensory Systems, Neuroscience), according to OpenAlex
Funding: National Institute on Deafness and Other Communication Disorders (R01DC016364, R01DC018539, R01DC021663); NINDS (T32NS047987)
Citations: cited by 1 paper (Europe PMC); 121 references in the paper
Research resources: RRID:SCR_001622, via the PsychToolBox extension RRID:SCR_002881, RRID:SCR_004849

Abstract

When sampling odors, humans typically take a single long sniff, unlike most other mammals who sample odors through bouts of rapid repetitive sniffing. Decades of work has established that rapid sniffing rhythms underlie the organizational principles of odor coding, with sniff speed clocking odor responses in the olfactory bulb. In the absence of rapid sniffing, how are odor responses organized in the human olfactory system? Because most mammals sniff at rates centered around the theta frequency, we hypothesized that the olfactory bulb exploits a theta-range neural oscillation to set the pace of odor coding in the system. Here, we used high-precision neural recordings from the human olfactory bulb to show that initiation of a sniff elicits and temporally aligns theta oscillations in the bulb and that sniff-induced theta oscillations organize the timing and amplitude of odor responses. These findings suggest that, despite the lack of rapid sniffing bouts in humans, the system has preserved a similarly timed unit of olfactory processing.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 16 matches between paragraphs and lines of code.

Zenodo 18791231

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Languages: MATLAB (22)
Size: 49 files, 22 scripts
Software Heritage: not checked
Found in: “Data, code, and materials availability:”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (7 files), FieldTrip (6 files), Image Processing Toolbox (6 files), Signal Processing Toolbox (5 files), shadedErrorBar (3 files), CircStat (2 files), Chronux (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
22 files
At the source: dx.doi.org/10.5281/zenodo.18791231

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 22 scripts, each with its path and the digest of its content;
  • 16 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data, code, and materials availability

All data and code needed to evaluate and reproduce the results in the paper are present in the paper and/or the Supplementary Materials. In addition, data and code are available on Zenodo (DOI: 10.5281/zenodo.18791231 (http://dx.doi.org/10.5281/zenodo.18791231), available at https://zenodo.org/records/18791231). This study did not generate new materials.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 13 authors, 6 MeSH terms, 2 funders, 118 references, 3 RRIDs.

Cite

This paper

Sheriff, A., Omidbeigi, M., Lane, G., Yang, Q., Zhou, G., Dede, A., Arriaga, N., Sagar, V., Holubecki, A. M., Braga, R. M., Kay, L. M., Tan, B. K., & Zelano, C. (2026). Theta oscillations are an organizational unit of odor processing in the olfactory bulb. Science advances, 12(27), eaee1002. https://doi.org/10.1126/sciadv.aee1002

BibTeX

@article{sheriff2026theta,
author = {Sheriff, Andrew and Omidbeigi, Mahmoud and Lane, Gregory and Yang, Qiaohan and Zhou, Guangyu and Dede, Adam and Arriaga, Naelly and Sagar, Vivek and Holubecki, Ania M. and Braga, Rodrigo M. and Kay, Leslie M. and Tan, Bruce K. and Zelano, Christina},
title = {{Theta oscillations are an organizational unit of odor processing in the olfactory bulb}},
journal = {Science advances},
year = {2026},
month = jul,
volume = {12},
number = {27},
pages = {eaee1002},
publisher = {American Association for the Advancement of Science},
issn = {2375-2548},
doi = {10.1126/sciadv.aee1002},
url = {https://doi.org/10.1126/sciadv.aee1002},
pmid = {42397924},
pmcid = {PMC13330860}
}

RIS

TY - JOUR
AU - Sheriff, Andrew
AU - Omidbeigi, Mahmoud
AU - Lane, Gregory
AU - Yang, Qiaohan
AU - Zhou, Guangyu
AU - Dede, Adam
AU - Arriaga, Naelly
AU - Sagar, Vivek
AU - Holubecki, Ania M.
AU - Braga, Rodrigo M.
AU - Kay, Leslie M.
AU - Tan, Bruce K.
AU - Zelano, Christina
TI - Theta oscillations are an organizational unit of odor processing in the olfactory bulb
T2 - Science advances
J2 - Sci Adv
PY - 2026
DA - 2026/07/03
VL - 12
IS - 27
SP - eaee1002
SN - 2375-2548
PB - American Association for the Advancement of Science
DO - 10.1126/sciadv.aee1002
UR - https://doi.org/10.1126/sciadv.aee1002
LA - en
ER -

CSL-JSON

{
"id": "10.1126/sciadv.aee1002",
"type": "article-journal",
"title": "Theta oscillations are an organizational unit of odor processing in the olfactory bulb",
"container-title": "Science advances",
"author": [
{
"family": "Sheriff",
"given": "Andrew"
},
{
"family": "Omidbeigi",
"given": "Mahmoud"
},
{
"family": "Lane",
"given": "Gregory"
},
{
"family": "Yang",
"given": "Qiaohan"
},
{
"family": "Zhou",
"given": "Guangyu"
},
{
"family": "Dede",
"given": "Adam"
},
{
"family": "Arriaga",
"given": "Naelly"
},
{
"family": "Sagar",
"given": "Vivek"
},
{
"family": "Holubecki",
"given": "Ania M."
},
{
"family": "Braga",
"given": "Rodrigo M."
},
{
"family": "Kay",
"given": "Leslie M."
},
{
"family": "Tan",
"given": "Bruce K."
},
{
"family": "Zelano",
"given": "Christina"
}
],
"container-title-short": "Sci Adv",
"volume": "12",
"issue": "27",
"page": "eaee1002",
"DOI": "10.1126/sciadv.aee1002",
"PMID": "42397924",
"PMCID": "PMC13330860",
"ISSN": "2375-2548",
"publisher": "American Association for the Advancement of Science",
"URL": "https://doi.org/10.1126/sciadv.aee1002",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
3
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-73106-z [code]
Respiratory pauses highlight sleep architecture in mice.
Journal: Nature communications
In common: Chronux, CircStat, FieldTrip, 3 other tools, 4 references
[2] doi:10.1371/journal.pbio.3003810 [code]
Olfactory bulb-cortex oscillations encode perceived odor intensity rather than concentration.
Journal: PLoS biology
In common: FieldTrip, Image Processing Toolbox, Statistics and Machine Learning Toolbox, 7 references
[3] doi:10.7554/elife.110795 [code]
REM sleep prefrontal high-frequency oscillation chains mediate distinct cortical - hippocampal reactivation patterns compared to NREM sleep.
Journal: eLife
In common: Chronux, shadedErrorBar, CircStat, 4 other tools, 3 references
[4] doi:10.1038/s41467-026-72445-1 [code]
Perception and neural representation of intermittent odor stimuli in mice.
Journal: Nature communications
In common: 8 references
[5] doi:10.64898/2026.03.12.710517 [code]
Cortical excitability inversely modulates fMRI connectivity via low-frequency neuronal coupling
Journal: bioRxiv (preprint)
In common: Chronux, shadedErrorBar, CircStat, 4 other tools, 1 reference
[6] doi:10.1038/s41467-026-73001-7 [code]
Separable and integrated pleasantness coding for appetitive and aversive odors across olfactory and ventral prefrontal cortices.
Journal: Nature communications
In common: Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 2 authors
[7] doi:10.1126/sciadv.aed3610 [code]
Dexterous single sniffs for ethological active olfaction.
Journal: Science advances
In common: 8 references
[8] doi:10.1016/j.isci.2026.116240 [code]
Physical activity enhances theta-periodicity of visual attentional allocation.
Journal: iScience
In common: shadedErrorBar, CircStat, Signal Processing Toolbox, 1 other tool, 4 references
[9] doi:10.1371/journal.pbio.3003938 [code]
Theta oscillations tag episodic memories for sleep-dependent consolidation.
Journal: PLoS biology
In common: shadedErrorBar, CircStat, FieldTrip, 3 other tools, 2 references
[10] doi:10.1016/j.celrep.2026.117646 [code]
Medial entorhinal-hippocampal desynchronization parallels the emergence of memory impairment in a mouse model of Alzheimer's disease pathology.
Journal: Cell reports
In common: Chronux, CircStat, FieldTrip, 3 other tools, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.