% Analyse dataset for closed-loop stimulation of EMX-ChR2 slices (Seizures) % Creates power spectra and modulation vs phase plots (Figure 6) % session() dataset contains 10 sessions % .name filename % .cl_freq closed-loop filter frequency (Hz) % .sample_rate sampling rate (Hz) % .num_sect number of data sections % .sect() datasection % .lfp local field potential (uV) % .stim stimulus (0-5) % .phase_cond phase condition % 1,2,..8 is 0,45..315 deg % 9 is no stimulation % .burst identified bursts % [duration, start time, phase condition, burst number] clear all load CLOSe_dataset_3 col=[0 0 1; 0 0.3 0.6; 0 0.6 0.3; 0 1 0; 0.3 0.6 0; 0.6 0.3 0; 1 0 0; 0.5 0 0.5; 0 0 0]; nfft=2^10; prange=[0:50]; srange=16:19; for dataset=1:length(session) freq=[0:nfft/2]*session(dataset).sample_rate/nfft; count=zeros(9,1); power_spec=zeros(nfft/2+1,9); band_power=zeros(9,1); band_power2=zeros(9,1); all_burst=[]; for i=1:session(dataset).num_sect % Compile power spectra temp=pwelch(session(dataset).sect(i).lfp,nfft,7*nfft/8,nfft,session(dataset).sample_rate); power_spec(:,session(dataset).sect(i).phase_cond)=power_spec(:,session(dataset).sect(i).phase_cond)+temp; % Calculate power modulation around closed-loop frequency temp=interp1(freq*20/session(dataset).cl_freq,temp,prange,'linear'); temp=sum(temp(srange)); band_power(session(dataset).sect(i).phase_cond)=band_power(session(dataset).sect(i).phase_cond)+temp; band_power2(session(dataset).sect(i).phase_cond)=band_power2(session(dataset).sect(i).phase_cond)+temp.^2; count(session(dataset).sect(i).phase_cond)=count(session(dataset).sect(i).phase_cond)+1; all_burst=[all_burst;session(dataset).sect(i).burst]; end % Plot power spectra figure subplot(2,2,1) hold on for c=9:-1:1 power_spec(:,c)=power_spec(:,c)/count(c); band_power(c)=band_power(c)/count(c); band_power2(c)=band_power2(c)/count(c); band_power2(c)=sqrt((band_power2(c)-band_power(c).^2)/(count(c)-1)); h=plot(freq,power_spec(:,c),'k'); if c<9 set(h,'Color',col(c,:)); end end xlabel('Frequency (Hz)') ylabel('Power (uV^2/Hz)') title([session(dataset).name]) axis([0 50 0 1.2*max(max(power_spec(5:50,:)))]) % Normalise power spectra power_modulation=[];band_power_modulation=[];band_power_modulation_se=[]; for c=1:8 band_power_modulation(c)=10*log10(band_power(c)./band_power(9)); band_power_modulation_se(c)=10*sqrt(band_power2(c).^2/band_power(c).^2+band_power2(9).^2/band_power(9).^2); power_modulation(:,c)=10*log10(power_spec(:,c)./power_spec(:,9)); end power_modulation(:,9:16)=power_modulation;power_modulation(:,17)=power_modulation(:,1); % Plot power modulation subplot(2,2,2) pcolor([-360:45:360],freq,power_modulation) colormap jet shading flat axis([-360 360 0 50]) caxis([-10 10]) xlabel('Phase (deg)') ylabel('Frequency (Hz)') set(gca,'XTick',[-360:180:360]); % Plot band-power modulation subplot(2,2,4) plot([-360:45:360],[band_power_modulation band_power_modulation band_power_modulation(1)],'b'); hold on plot([-360:45:360],[band_power_modulation band_power_modulation band_power_modulation(1)]+[band_power_modulation_se band_power_modulation_se band_power_modulation_se(1)],'b--'); plot([-360:45:360],[band_power_modulation band_power_modulation band_power_modulation(1)]-[band_power_modulation_se band_power_modulation_se band_power_modulation_se(1)],'b--'); set(gca,'XLim',[-360 360]) set(gca,'XTick',[-360:180:360]); xlabel('Phase (deg)') ylabel('Power modulation (dB)') % Calculate relative burst duration (separately for first and subsequent bursts) first_burst=all_burst(find(all_burst(:,4)==1),[1 3]); subsequent_burst=all_burst(find(all_burst(:,4)~=1),[1 3]); no_stim=mean(log2(first_burst(find(first_burst(:,2)==9),1))); no_stim(2)=mean(log2(subsequent_burst(find(subsequent_burst(:,2)==9),1))); first_burst(:,1)=log2(first_burst(:,1))-no_stim(1); subsequent_burst(:,1)=log2(subsequent_burst(:,1))-no_stim(2); all_burst=[first_burst;subsequent_burst]; burst_duration_modulation=[];burst_duration_modulation_se=[]; for i=1:8 burst_duration_modulation(i)=mean(all_burst(find(all_burst(:,2)==i),1)); burst_duration_modulation_se(i)=std(all_burst(find(all_burst(:,2)==i),1))/sqrt(length(find(all_burst(:,2)==i))); end burst_duration_modulation=burst_duration_modulation';burst_duration_modulation_se=burst_duration_modulation_se'; % Plot relative burst duration subplot(2,2,3) plot([-8:8]*45,2.^[burst_duration_modulation;burst_duration_modulation;burst_duration_modulation(1)],'b') hold on plot([-8:8]*45,2.^[burst_duration_modulation;burst_duration_modulation;burst_duration_modulation(1)]+[burst_duration_modulation_se;burst_duration_modulation_se;burst_duration_modulation_se(1)],'b--') plot([-8:8]*45,2.^[burst_duration_modulation;burst_duration_modulation;burst_duration_modulation(1)]-[burst_duration_modulation_se;burst_duration_modulation_se;burst_duration_modulation_se(1)],'b--') xlabel('Phase (deg)') ylabel('Relative burst duration') set(gca,'XTick',[-360:180:360]); set(gca,'XLim',[-360 360]); drawnow end