%%%% script to display graphs from models-as-animals data (LFO single cell fitting)
clear all
%%
exp_no = 18;
freq_display_limit = 40;
%% paths
batch_path = '../ResultsArchive/LFO-urethane/CompleteModel/ConditionA/';
%batch_path = '../ResultsArchive/LFO-urethane/CompleteModel/ConditionB/';
%batch_path = '../ResultsArchive/LFO-urethane/CompleteModel/ConditionC/';
%batch_path = '../ResultsArchive/LFO-urethane/CompleteModel/ConditionD/';
% Collaterals, no STN 3rd DA: batch of 50 - Normal model
batch = 'LFO_a_20060407T004326_batch.mat';
%batch = 'LFO_b_20060407T011449_batch.mat';
%batch = 'LFO_c_20060407T003124_batch.mat';
%batch = 'LFO_d_20060407T000419_batch.mat';
% parameters for analyses done here
ISIbins = 100; %No. of bins in ISI hist.
limits = [0.005 0.6]; % ISI hist limits in seconds (outside these bounds get put into end bins)
trig_win_half = 1;
acf_binsize = 0.01; % in seconds
acf_maxlag = 2; % window size for auto-correlograms
% parameters for periodogram
freq_range = [0.1 10]; % min and max freqs in periodogram
tapers = [3 5]; % recommended in Chronuz; 5 tapers with bandwidth*time = 3
pad = 2; % Padding the frequency sampling by another 2 powers of 2;
err_bars = [1 0.01]; % method 1 (theoretical) and p_sig = 0.01
fscorr = 0; %don't use finite size corrections
%% set paths for interactive sessions
[a host] = system('hostname');
%%% set requisite paths!
if findstr(host, 'iceberg'); % on iceberg
fprintf('\n On ICEBERG \n');
system_os = 'unix';
ice_path1 = genpath('/home1/pc/pc1mdh/BG spiking model');
ice_path2 = genpath('/home1/pc/pc1mdh/Matlab Tools');
path(path, ice_path1);
path(path, ice_path2);
elseif (findstr(host, 'node') | findstr(host,'ace')) % on ACE
system_os = 'unix';
ace_path1 = genpath('/home/mark/SpikingModel');
path(path, ace_path1);
fprintf('\n On ACE \n');
else
system_os = 'xp';
fprintf('\n On XP \n');
end
%% load batch lists find the number of batches
load([batch_path batch]);
[n_batches c] = size(batch_analysis_list);
%%%%%%%%%% do Condition A (control) %%%%%%%%%%%%%%
extracted_list = batch_analysis_list{exp_no,3};
analysis_list= batch_analysis_list{exp_no,2};
n_animals = length(analysis_list);
cl_fig
for loop1 = 1:n_animals
load([batch_path extracted_list{loop1}])
load([batch_path analysis_list{loop1}])
fprintf(1, 'mean STN rate %.2f: stderr of mean STN rate %.2f: std dev of STN rate %.2f\n', STN_Hz, sem_STN, std_STN);
fprintf(1, 'mean GPe rate %.2f: stderr of mean GPe rate %.2f: std dev of STN rate %.2f\n', GPe_Hz, sem_GPe, std_GPe);
fprintf(1, 'mean GPi rate %.2f: stderr of mean GPi rate %.2f: std dev of STN rate %.2f\n', GPi_Hz, sem_GPi, std_GPi);
% do cortical EEG power spectrum
new_dt = dt * 10;
time_stamps = new_dt:new_dt:time_seconds;
[ctx_result,ctx_power,ctx_fs] = scargle_analysis(ctx_mean,time_stamps,new_dt,time_seconds,256,[0.1 5]);
%%%%%%%%%% show data for STN %%%%%%%%%
STN_acf = cell(n_STN,1);
xSTNacf = cell(n_STN,1);
temp_acf = cell(n_STN,1);
limit = cell(n_STN,1);
% do power spectrum of cortical wave
% Fs = 1 ./ dt; % for spectrum best to use the underlying sampling;
% [ctx_MT_power, ctx_MT_fs_array, R, ctx_MT_errs] = mtspectrumpt(ctx_mean', tapers, pad, Fs, freq_range, err_bars, 0, fscorr);
for loop = 1:n_STN
figure
subplot(2,3,1), plot(ctx_mean);
[STN_rates{loop},STN_ISIhist{loop},STN_x_isi{loop},STN_x_hist] = LIF_ISI_analysis(STN_times{loop},ISIbins,limits);
subplot(2, 3, 4);
plot(STN_x_isi{loop}, STN_rates{loop});
subplot(2,3,2), plot(ctx_fs,ctx_power)
% run auto-corr, take out centre bin
[STN_acf{loop},xSTNacf{loop},f1,f2,n_pairs] = LIF_xcorr(STN_times{loop},STN_times{loop}, acf_binsize, [0 time_seconds], acf_maxlag);
temp_acf{loop} = STN_acf{loop};
temp_acf{loop}(find(xSTNacf{loop}>=-0.05 & xSTNacf{loop}<= 0.05)) = mean(STN_acf{loop});
subplot(2,3,5), plot(xSTNacf{loop},temp_acf{loop})
% plot cell spectra 0-10Hz
limit{loop} = find(STN_spect_res(loop).freqs <= freq_display_limit);
subplot(2,3,3), plot(STN_spect_res(loop).freqs(limit{loop}),STN_spect_res(loop).powers(limit{loop}));
% make time-scale for x-axis (+/- 1 s?)
trig_x = linspace(-trig_win_half,trig_win_half,length(STN_avg_wvfrm(loop,:)));
% plot spike-triggered average waveform
subplot(2,3,6), plot(trig_x,STN_avg_wvfrm(loop,:));
h = gca;
txt = ['animal ' num2str(loop1) '- STN, cell ' num2str(loop)];
title(txt);
drawnow
end
tile
resp = input('save STN cell data for? (enter each number e.g. 41, enter 0 to skip)');
if resp == 0
else
resp = num2str(resp);
for i = 1:length(resp);
cell_num = str2num(resp(i));
file_prefix = [batch(1:end-9) '_exp_' num2str(exp_no) '_animal_' num2str(loop1) '_STNcell_' resp(i)];
% ctx eeg
ctx_eeg = [time_stamps' ctx_mean'];
save([file_prefix '_ctx_eeg.txt'],'ctx_eeg','-ascii');
% generate spikes before saving
ts = STN_times{cell_num};
marks = -60 .* ones(length(ts), 1);
stn_ras = [ts' marks];
save([file_prefix '_ras.txt'],'stn_ras','-ascii');
% cortex Lomb
ctx_Lomb = [ctx_fs' ctx_power'];
save([file_prefix '_ctx_Lomb.txt'],'ctx_Lomb','-ascii');
% autocorrelogram
stn_acorr = [xSTNacf{cell_num}' temp_acf{cell_num}'];
save([file_prefix '_acorr.txt'],'stn_acorr','-ascii');
% power spectrum
stn_power = [STN_spect_res(cell_num).freqs' STN_spect_res(cell_num).powers];
save([file_prefix '_power.txt'],'stn_power','-ascii');
% spike-triggered avg waveform
stn_wvfrm = [trig_x' STN_avg_wvfrm(cell_num,:)'];
save([file_prefix '_avg_wvfrm.txt'],'stn_wvfrm','-ascii');
end
end
cl_fig
% and the same for the GPe.....
GPe_acf = cell(n_GPe,1);
xGPeacf = cell(n_GPe,1);
temp_acf = cell(n_GPe,1);
limit = cell(n_GPe,1);
for loop = 1:n_GPe
figure
subplot(2,3,1), plot(ctx_mean);
[GPe_rates{loop},GPe_ISIhist{loop},GPe_x_isi{loop},GPe_x_hist] = LIF_ISI_analysis(GPe_times{loop},ISIbins,limits);
subplot(2, 3, 4);
plot(GPe_x_isi{loop}, GPe_rates{loop});
subplot(2,3,2), plot(ctx_fs,ctx_power)
% run auto-corr, take out centre bin
[GPe_acf{loop},xGPeacf{loop},f1,f2,n_pairs] = LIF_xcorr(GPe_times{loop},GPe_times{loop}, acf_binsize, [0 time_seconds], acf_maxlag);
temp_acf{loop} = GPe_acf{loop};
temp_acf{loop}(find(xGPeacf{loop}>=-0.05 & xGPeacf{loop} <= 0.05)) = mean(GPe_acf{loop});
subplot(2,3,5), plot(xGPeacf{loop},temp_acf{loop})
% plot cell spectra 0-10Hz
limit{loop} = find(GPe_spect_res(loop).freqs <= freq_display_limit);
subplot(2,3,3), plot(GPe_spect_res(loop).freqs(limit{loop}),GPe_spect_res(loop).powers(limit{loop}));
% make time-scale for x-axis (+/- 1 s?)
trig_x = linspace(-trig_win_half,trig_win_half,length(GPe_avg_wvfrm(loop,:)));
% plot spike-triggered average waveform
subplot(2,3,6), plot(trig_x,GPe_avg_wvfrm(loop,:));
h = gca;
txt = ['animal ' num2str(loop1) '- GPe, cell ' num2str(loop)];
title(txt);
drawnow
end
tile
resp = input('save GPe cell data for? (enter each number e.g. 41, enter 0 to skip)');
if resp == 0
else
resp = num2str(resp);
for i = 1:length(resp);
cell_num = str2num(resp(i));
file_prefix = [batch(1:end-9) '_exp_' num2str(exp_no) '_animal_' num2str(loop1) '_GPecell_' resp(i)];
% ctx eeg
ctx_eeg = [time_stamps' ctx_mean'];
save([file_prefix '_ctx_eeg.txt'],'ctx_eeg','-ascii');
% generate spikes before saving
ts = GPe_times{cell_num};
marks = -60 .* ones(length(ts), 1);
stn_ras = [ts' marks];
save([file_prefix '_ras.txt'],'stn_ras','-ascii');
% cortex Lomb
ctx_Lomb = [ctx_fs' ctx_power'];
save([file_prefix '_ctx_Lomb.txt'],'ctx_Lomb','-ascii');
% autocorrelogram
stn_acorr = [xGPeacf{cell_num}' temp_acf{cell_num}'];
save([file_prefix '_acorr.txt'],'stn_acorr','-ascii');
% power spectrum
stn_power = [GPe_spect_res(cell_num).freqs' GPe_spect_res(cell_num).powers];
save([file_prefix '_power.txt'],'stn_power','-ascii');
% spike-triggered avg waveform
stn_wvfrm = [trig_x' GPe_avg_wvfrm(cell_num,:)'];
save([file_prefix '_avg_wvfrm.txt'],'stn_wvfrm','-ascii');
end
end
cl_fig
end