% Plotting the hemodynamic response function (hrf) using fmridesign hrf_parameters=[5.4 5.2 10.8 7.35 0.35 0] time=(0:240)/10; hrf0=fmridesign(time,0,[1 0],[],hrf_parameters); plot(time,hrf0,'LineWidth',2) xlabel('time (seconds)') ylabel('hrf') title('Glover hrf model') % Making the design matrices using fmridesign frametimes=(0:119)*3; slicetimes=[0.14 0.98 0.26 1.10 0.38 1.22 0.50 1.34 0.62 1.46 0.74 1.58 0.86]; eventid=kron(ones(10,1),[1; 2]); eventimes=(0:19)'*18+9; duration=ones(20,1)*9; height=ones(20,1); events=[eventid eventimes duration height] hrf_parameters=[5.4 5.2 10.8 7.35 0.35 0] S=kron(ones(10,1),kron([0 0; 1 0; 0 0; 0 1],ones(3,1))); X_cache=fmridesign(frametimes,slicetimes,events,[],hrf_parameters); X_cache=fmridesign(frametimes,slicetimes, [] , S,hrf_parameters); plot(X_cache(:,9:10),'LineWidth',2) legend('Hot','Warm') xlabel('frame number') ylabel('response') title('Hat and Warm responses') % Analysing one run with fmrilm contrast=[1 0 0 0 0 0; 0 1 0 0 0 0; 1 -1 0 0 0 0]; exclude=[1 2]; which_stats=[1 1 1 1]; contrast=[0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1] % A linear effect of temperature temperature=[45.5 35.0 49.0 38.5 42.0 49.0 35.0 42.0 38.5 45.5 ... 38.5 49.0 35.0 45.5 42.0 45.5 38.5 42.0 35.0 49.0]'; Temperature=35:3.5:49; subplot(2,2,1) plot(Temperature,1:5,'LineWidth',2); xlabel('Temperature'); ylabel('Response'); title('(a) Linear temperature effect'); subplot(2,2,2) plot(Temperature,10-((1:5)-4).^2,'LineWidth',2); xlabel('Temperature'); ylabel('Response'); title('(b) Quadratic temperature effect'); subplot(2,2,3) plot(Temperature,10-((1:5)-4).^2+((1:5)-3).^3,'LineWidth',2); xlabel('Temperature'); ylabel('Response'); title('(c) Cubic temperature effect'); subplot(2,2,4) plot(Temperature,[1 4 3 5 3.5],'LineWidth',2); xlabel('Temperature'); ylabel('Response'); title(['(d) Quartic = arbitrary temperature effect']); events=[zeros(20,1)+1 eventimes duration ones(20,1); zeros(20,1)+2 eventimes duration temperature] contrast=[0 1 0 0 0 0; 1 49 0 0 0 0; 1 35 0 0 0 0; 0 14 0 0 0 0]; events=[zeros(20,1)+1 eventimes duration ones(20,1); zeros(20,1)+2 eventimes duration temperature; zeros(20,1)+3 eventimes duration temperature^2] contrast=[0 0 1 0 0 0 0]; contrast=[1 0 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0 0; 0 0 1 0 0 0 0 0 0; 0 0 0 1 0 0 0 0 0; 0 0 0 0 1 0 0 0 0]; contrast=[eye(5) ones(5,4)]; which_stats=[0 0 0 1]; contrast=[.8 -.2 -.2 -.2 -.2 0 0 0 0; -.2 .8 -.2 -.2 -.2 0 0 0 0; -.2 -.2 .8 -.2 -.2 0 0 0 0; -.2 -.2 -.2 .8 -.2 0 0 0 0; -.2 -.2 -.2 -.2 .8 0 0 0 0]; contrast=[eye(5)-ones(5)/5 ones(5,4)]; which_stats=[0 0 0 1]; contrast=[-7.0 -3.5 0 3.5 7.0 0 0 0 0]; which_stats=[1 1 1]; % Combining runs/sessions/subjects with multistat input_file='c:\keith\data\brian_ha_19971124_1_093923_mri_MC.mnc'; output_file_base=['c:\keith\data\ha_093923_hot'; 'c:\keith\data\ha_093923_wrm'; 'c:\keith\data\ha_093923_hmw']; [df1 p]=fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats); input_file='c:\keith\data\brian_ha_19971124_1_100326_mri_MC.mnc'; output_file_base=['c:\keith\data\ha_100326_hot'; 'c:\keith\data\ha_100326_wrm'; 'c:\keith\data\ha_100326_hmw']; [df2 p]=fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats); input_file='c:\keith\data\brian_ha_19971124_1_101410_mri_MC.mnc'; output_file_base=['c:\keith\data\ha_101410_hot'; 'c:\keith\data\ha_101410_wrm'; 'c:\keith\data\ha_101410_hmw']; [df3 p]=fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats); input_file='c:\keith\data\brian_ha_19971124_1_102703_mri_MC.mnc'; output_file_base=['c:\keith\data\ha_102703_hot'; 'c:\keith\data\ha_102703_wrm'; 'c:\keith\data\ha_102703_hmw']; [df4 p]=fmrilm(input_file, output_file_base, X_cache, contrast, exclude, which_stats); X=[1 1 1 1]'; contrast=[1]; which_stats=[1 1 1]; input_files_effect=['c:\keith\data\ha_093923_hmw_effect.mnc'; 'c:\keith\data\ha_100326_hmw_effect.mnc'; 'c:\keith\data\ha_101410_hmw_effect.mnc'; 'c:\keith\data\ha_102703_hmw_effect.mnc']; input_files_sdeffect=['c:\keith\data\ha_093923_hmw_sdeffect.mnc'; 'c:\keith\data\ha_100326_hmw_sdeffect.mnc'; 'c:\keith\data\ha_101410_hmw_sdeffect.mnc'; 'c:\keith\data\ha_102703_hmw_sdeffect.mnc']; output_file_base='c:\keith\data\ha_multi_hmw' input_files_df=[df1 df2 df3 df4] input_files_fwhm=6 df=multistat(input_files_effect,input_files_sdeffect,input_files_df, ... input_files_fwhm,X,contrast,output_file_base,which_stats,Inf) % Fixed and random effects df=multistat(input_files_effect,input_files_sdeffect,input_files_df, ... input_files_fwhm,X,contrast,output_file_base,which_stats,15) % Thresholding the tstat image with tstat_threshold tstat_threshold(1000000,38.4521,6,112,0.05) tstat_threshold(1000000,38.4521,6,112,0.05,3) % Locating peaks with locmax lm=locmax([output_file_base '_tstat.mnc'],thresh) pval=tstat_threshold(1000000,38.4521,6,112,lm(:,1)) % Extracting values from a minc file using extract voxel=[62 67 4] effect=extract(voxel,output_file_base,'_effect.mnc') sdeffect=extract(voxel,output_file_base,'_sdeffect.mnc') % Estimating the delay hrf_parameters2=[5.4 5.2 10.8 7.35 0.35 2] X_cache2=fmridesign(frametimes,slicetimes,events,[],hrf_parameters2); contrast=[1 0 0 0 0 0; 0 1 0 0 0 0] which_stats=[1 1 1 0] input_file='/scratch/keith/brian_ha_19971124_1_100326_mri_MC.mnc'; output_file_base=['/scratch/keith/ha_100326_delay_hot'; '/scratch/keith/ha_100326_delay_wrm'] n_poly=3; fwhm_rho=15; reference_delay=5.4; fmrilm(input_file, output_file_base, X_cache2, contrast, ... exclude, which_stats, fwhm_rho, n_poly, reference_delay); % Estimating the time course of the response eventid=kron(ones(10,1),(1:12)'); eventimes=frametimes'; duration=ones(120,1)*3; height=ones(120,1); events=[eventid eventimes duration height] X_cache=fmridesign(frametimes,slicetimes,events,[],zeros(1,6)); plot(X_cache(:,4*10+(1:4)),'LineWidth',2) legend('Time01','Time02','Time03','Time04') xlabel('frame number') ylabel('response') title('Slice 5, first 4 times only') contrast=[eye(12)-ones(12)/12 zeros(12,4)]; exclude=[1 2]; which_stats=[0 1 1 1]; input_file='c:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc'; output_file_base=['c:/keith/data/ha_100326_time01'; 'c:/keith/data/ha_100326_time02'; 'c:/keith/data/ha_100326_time03'; 'c:/keith/data/ha_100326_time04'; 'c:/keith/data/ha_100326_time05'; 'c:/keith/data/ha_100326_time06'; 'c:/keith/data/ha_100326_time07'; 'c:/keith/data/ha_100326_time08'; 'c:/keith/data/ha_100326_time09'; 'c:/keith/data/ha_100326_time10'; 'c:/keith/data/ha_100326_time11'; 'c:/keith/data/ha_100326_time12']; fmrilm(input_file,output_file_base,X_cache,contrast,exclude,which_stats); tstat_threshold(1000000,38.4538,6,[11 103]) lm=locmax([output_file_base(1,:) '_fstat.mnc'],5.18) voxel=[63 66 4] values=extract(voxel,output_file_base,'_effect.mnc') sd=extract(voxel,output_file_base,'_sdeffect.mnc') b_hot=extract(voxel,'c:/keith/data/ha_100326_hot','_effect.mnc') b_wrm=extract(voxel,'c:/keith/data/ha_100326_wrm','_effect.mnc') time=(1:360)/10; hrf_parameters=[5.4 5.2 10.8 7.35 0.35 0] hrf=fmridesign(time,0,[1 9 9 1],[],hrf_parameters); plot((0:12)*3+slicetimes(voxel(3)),values([1:12 1]),'k', ... [0:11; 0:11]*3+slicetimes(voxel(3)), [values+sd; values-sd],'g', ... time,[zeros(1,90) ones(1,90) zeros(1,180)],'r', ... time,[zeros(1,270) ones(1,90)],'b', ... time,hrf*b_hot+hrf([181:360 1:180])*b_wrm,'g','LineWidth',2); legend('Estimated response','Modeled response'); text(10,0.5,'Hot') text(28,0.5,'Warm') xlabel('time (seconds) from start of epoch'); ylabel('fMRI response'); % Estimating the hrf hrf_peaks=[3:3:12] np=length(hrf_peaks) hrf_parameters=[ [hrf_peaks hrf_peaks]' 3*ones(2*np,1) zeros(2*np,4)] time=(0:240)/10; hrfs=fmridesign(time,0,[(1:np)' zeros(np,1)],[],hrf_parameters(1:np,:)); plot(time,hrfs,'LineWidth',2); xlabel('time (seconds)'); ylabel('hrf'); title('Hrf basis functions'); eventid=kron(ones(10,1),(1:2*np)'); eventimes=kron((0:19)'*18+9,ones(np,1)); duration=ones(20*np,1)*9; height=ones(20*np,1); events=[eventid eventimes duration height] X_cache=fmridesign(frametimes,slicetimes,events,[],hrf_parameters); plot(X_cache(1:24,1:2*np),'LineWidth',2) xlabel('frame number'); ylabel('fMRI response'); title('Hrf basis functions convolved with the hot and warm stimuli') contrast=[eye(2*np)-ones(2*np)/2/np zeros(2*np,4)] which_stats=[0 1 1 1]; input_file='c:/keith/data/brian_ha_19971124_1_100326_mri_MC.mnc'; output_file_base=['c:/keith/data/ha_100326_hot1'; 'c:/keith/data/ha_100326_hot2'; 'c:/keith/data/ha_100326_hot3'; 'c:/keith/data/ha_100326_hot4'; 'c:/keith/data/ha_100326_wrm1'; 'c:/keith/data/ha_100326_wrm2'; 'c:/keith/data/ha_100326_wrm3'; 'c:/keith/data/ha_100326_wrm4']; fmrilm(input_file, output_file_base, X_cache,contrast,exclude, which_stats); tstat_threshold(1000000,38.4538,6,[7 106]) lm=locmax([output_file_base(1,:) '_fstat.mnc'],6.5642) voxel=[62 65 4] coefs=extract(voxel,output_file_base,'_effect.mnc') sdcoefs=extract(voxel,output_file_base,'_sdeffect.mnc') b_hot=extract(voxel,'c:/keith/data/ha_100326_hot','_effect.mnc') b_wrm=extract(voxel,'c:/keith/data/ha_100326_wrm','_effect.mnc') hrf0=fmridesign(time); plot(time, [hrfs*coefs(1:np)' hrf0*b_hot hrfs*coefs((1:np)+np)' hrf0*b_wrm], ... 'LineWidth',2); legend('Hot estimated','Hot model','Warm estimated','Warm model'); xlabel('time (seconds)'); ylabel('hrf'); title('Estimated hrfs for the hot and warm stimuli, together with the Glover hrfs')