%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RESULT ANALYSIS FOR EBCC SIMULATION ISI = 300 ms %
% ------------------- %
% copyright : (C) 2014 by Alberto Antonietti %
% email : alberto.antonietti@polimi.it %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INSTRUCTIONS %
% ---------------- %
% This MATLAB file can be used to analyze the results of a simulation %
% of the Eye Blinking Classical Conditioning protocol using a SNN in %
% the EDLUT platform. Simply run this file after the simulation end %
% it will load the output file, the weight's files and the CR file. %
% This script generates different plots: %
% %
% - Raster Plots of the spikes generated by the network in 3 %
% different phases of the simulation (Beginning: 1st trial, %
% CR Acquired: 380th trial, Extinction: 580th trial - Fig. %
% 5 A-B-C left panels) %
% - Firing Rate Plots of the various cells during the same phases %
% (Fig. 5 A-B-C right panels) %
% - 3D Firing Rate Plots for PC and DCN for the whole simulation %
% (Fig. 6 A-B) %
% - CR Percentage computed with a mobile window of 10 trials %
% (Fig. 6 D) %
% - Synaptic Weights (PF-PC) Trend during the whole simulation %
% %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% RASTER PLOTS
clc
clear
load 'output.dat';
load 'DCNAvg.dat';
load('CR.dat');
repduration=500;
nrep=numel(DCNAvg)/repduration;
spk=sort(output(:,2)');
nspk=[];
if length(spk) > 0
reps=1;
nspk(1)=spk(1);
else
reps=0;
end
for n=2:length(spk),
if spk(n-1) ~= spk(n)
reps=reps+1;
nspk(reps)=spk(n);
end
end
CRs=CR;
for i=2:numel(CRs)
if CRs(i)<(CRs(i-1)+0.500)
CRs(i)=NaN;
end
end
CRs(find(isnan(CRs)))=[];
CRs=ceil(CRs.*1000);
counterCSs=1;
counterUSs=1;
counterCRs=1;
CR_all=ones(1,numel(DCNAvg))*NaN;
for i=1:size(CRs,1)
CR_all(CRs(counterCRs))=1;
counterCRs=counterCRs+1;
end
CR_all(find(CR_all==0))=NaN;
maxDCN=max(DCNAvg);
minDCN=min(DCNAvg);
if maxDCN==0
maxDCN=1;
end
for i = [0.001 190.5 290]
figure
set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.
starting = i;
ending = i+0.5;
h1=subplot(2,1,1);
if size(starting:0.001:ending)==size(starting*1000:ending*1000)
plot(starting:0.001:ending,DCNAvg(starting*1000:ending*1000),'k','LineWidth',2)
plot(starting:0.001:ending,DCNAvg(starting*1000:ending*1000),'k','LineWidth',2)
hold on
plot(starting:0.001:ending,CR_all(starting*1000:ending*1000).*maxDCN,'.k','Marker','pentagram','MarkerFaceColor',[0 0 0],'MarkerSize',10);
else
plot(starting:0.001:ending-0.001,DCNAvg(starting*1000:ending*1000),'k','LineWidth',2)
hold on
plot(starting:0.001:ending-0.001,CR_all(starting*1000:ending*1000).*maxDCN,'.k','Marker','pentagram','MarkerFaceColor',[0 0 0],'MarkerSize',10);
end
ylabel('Mean DCN Firing Rate [Hz]','FontSize',13);
ax=[0.0555555555555556 0.734579439252338 0.8986111111111111 0.244626168224301];
line([starting ending],[20 20],'Color',[1 0 0],'LineWidth',1.5);
ylim([minDCN*0.9 maxDCN*1.1])
xlim([starting ending])
box off
set(gca,'XTick',[]);
set(gca,'FontSize',15)
set(h1,'Position',ax);
h2=subplot(2,1,2);
ax=[0.0555555555555556 0.0355140186915888 0.8986111111111111 0.6959674627898926];
set(h2,'Position',ax);
output2=output(find(output(:,1)>=starting),:);
output2=output2(find(output2(:,1)<=ending),:);
tot_spks=0;
for n=1:reps
tspk=output2(find(output2(:,2)==nspk(n)),1);
tot_spks=tot_spks+length(tspk);
%MF
if(nspk(n)<=19)
line((tspk*ones(1,2))',(ones(length(tspk),1)*[n-0.43,n+0.43])','Color','b','LineWidth',2);
%PC
elseif(nspk(n)>=1520 && nspk(n)<=1543)
line((tspk*ones(1,2))',(ones(length(tspk),1)*[n-0.43,n+0.43])','Color','g','LineWidth',2);
%DCN
elseif(nspk(n)>=1544 && nspk(n)<=1555)
line((tspk*ones(1,2))',(ones(length(tspk),1)*[n-0.43,n+0.43])','Color','k','LineWidth',2);
%IO
elseif(nspk(n)>=1559 && nspk(n)<=1582)
line((tspk*ones(1,2))',(ones(length(tspk),1)*[n-0.43,n+0.43])','Color','m','LineWidth',2);
end
end
axis tight
xlabel('Time [s]');
display(['Total number of spikes: ' num2str(tot_spks)]);
display(['Number of spiking neurons: ' num2str(reps)]);
set(gca,'YTick',[]);
axis([starting ending 0 reps])
annotation('textbox',...
[0.0121111111111111 0.606024096385542 0.0253888888888889 0.048192771084338],...
'String',{'IO'},...
'FontSize',14,...
'FitBoxToText','off',...
'EdgeColor','none',...
'Color',[1 0 1]);
annotation('textbox',...
[0.0121111111111111 0.449830066534796 0.0253888888888889 0.048192771084338],...
'String',{'DCN'},...
'FontSize',14,...
'FitBoxToText','off',...
'EdgeColor','none');
annotation('textbox',...
[0.0121111111111111 0.307516633698974 0.0253888888888889 0.048192771084338],...
'String',{'PC'},...
'FontSize',14,...
'FitBoxToText','off',...
'EdgeColor','none',...
'Color',[0 1 0]);
annotation('textbox',...
[0.0121111111111111 0.111621111310914 0.0253888888888889 0.048192771084338],...
'String',{'MF'},...
'FontSize',14,...
'FitBoxToText','off',...
'EdgeColor','none',...
'Color',[0 0 1]);
set(gca,'FontSize',15)
annotation('textbox',...
[0.0081 0.005621111310914 0.0253888888888889 0.048192771084338],...
'String',{'Time'},...
'FontSize',14,...
'FitBoxToText','off',...
'EdgeColor','none',...
'Color',[0 0 0],'HorizontalAlignment','center');
set(gca,'FontSize',15)
end
%% FIRING RATE PLOTS
FirstMF=0;MFNum=20;
FirstGR=20;GRNum=1500;
FirstPC=1520;PCNum=24;
FirstDCN=1544;DCNNum=12;
FirstIO=1559;IONum=24;
window=25;
mobile=50;
maxMF=60;
maxPC=250;
maxDCN=30;
maxIO=20;
spk=sort(output(:,2)');
nspk=[];
if length(spk) > 0
reps=1;
nspk(1)=spk(1);
else
reps=0;
end
for n=2:length(spk),
if spk(n-1) ~= spk(n)
reps=reps+1;
nspk(reps)=spk(n);
end
end
for starting=[0.001 190.5 290]
ending = starting+0.5;
output_cut=output(output(:,1)>=starting,:);
output_cut=output_cut(output_cut(:,1)<ending,:);
% MF FREQUENCY ANALYSIS
fspk=zeros((ending-starting)*1000,MFNum);
for i=1:MFNum
for j=1:size(fspk,1)
tspk=output_cut(output_cut(:,2)==i-1+FirstMF,1);
tspk=ceil(tspk.*1000);
inizio=j;
fine=j+window;
a=tspk(tspk>=inizio+starting*1000);
a=a(a<fine+starting*1000);
fspk(j,i)=numel(a)*1000/window;
clear a
end
end
MF_meanfspk=smooth(mean(fspk,2),mobile);
clear fspk
% PC FREQUENCY ANALYSIS
fspk=zeros((ending-starting)*1000,PCNum);
for i=1:PCNum
for j=1:size(fspk,1)
tspk=output_cut(output_cut(:,2)==i-1+FirstPC,1);
tspk=ceil(tspk.*1000);
inizio=j;
fine=j+window;
a=tspk(tspk>=inizio+starting*1000);
a=a(a<fine+starting*1000);
fspk(j,i)=numel(a)*1000/window;
clear a
end
end
PC_meanfspk1=smooth(mean(fspk(:,1:PCNum/2),2),mobile);
PC_meanfspk2=smooth(mean(fspk(:,PCNum/2+1:PCNum),2),mobile);
clear fspk
% IO FREQUENCY ANALYSIS
fspk=zeros((ending-starting)*1000,IONum);
for i=1:IONum
for j=1:size(fspk,1)
tspk=output_cut(output_cut(:,2)==i-1+FirstIO,1);
tspk=ceil(tspk.*1000);
inizio=j;
fine=j+window;
a=tspk(tspk>=inizio+starting*1000);
a=a(a<fine+starting*1000);
fspk(j,i)=numel(a)*1000/window;
clear a
end
end
IO_meanfspk1=smooth(mean(fspk(:,1:IONum/2),2),mobile);
IO_meanfspk2=smooth(mean(fspk(:,IONum/2+1:IONum),2),mobile);
clear fspk
% DCN FREQUENCY ANALYSIS
fspk=zeros((ending-starting)*1000,DCNNum);
for i=1:DCNNum
for j=1:size(fspk,1)
tspk=output_cut(output_cut(:,2)==i-1+FirstDCN,1);
tspk=ceil(tspk.*1000);
inizio=j;
fine=j+window;
a=tspk(tspk>=inizio+starting*1000);
a=a(a<fine+starting*1000);
fspk(j,i)=numel(a)*1000/window;
clear a
end
end
DCN_meanfspk1=smooth(mean(fspk(:,1:DCNNum/2),2),mobile);
DCN_meanfspk2=smooth(mean(fspk(:,DCNNum/2+1:DCNNum),2),mobile);
clear fspk
% PLOT
for startingzoom = starting %SET THE START OF THE TRIAL
endingzoom=ending; %SET THE PERIOD OF THE TRIAL
% SUBPLOT
figure
set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.
h1=subplot(4,1,1);
plot(startingzoom:0.001:endingzoom-0.001,IO_meanfspk1,'m','LineWidth',3)
box off
hold on
plot(startingzoom:0.001:endingzoom-0.001,IO_meanfspk2,'--m','LineWidth',3)
set(gca,'XTick',[],'FontSize',16);
ylim([0 maxIO])
xlim([startingzoom endingzoom])
h2=subplot(4,1,2);
plot(startingzoom:0.001:endingzoom-0.001,DCN_meanfspk1,'k','LineWidth',3)
box off
hold on
plot(startingzoom:0.001:endingzoom-0.001,DCN_meanfspk2,'--k','LineWidth',3)
set(gca,'XTick',[],'FontSize',16);
ylim([0 maxDCN])
xlim([startingzoom endingzoom])
h3=subplot(4,1,3);
plot(startingzoom:0.001:endingzoom-0.001,PC_meanfspk1,'g','LineWidth',3)
box off
hold on
plot(startingzoom:0.001:endingzoom-0.001,PC_meanfspk2,'--g','LineWidth',3)
set(gca,'XTick',[],'FontSize',16);
ylim([0 maxPC])
xlim([startingzoom endingzoom])
h4=subplot(4,1,4);
plot(startingzoom:0.001:endingzoom-0.001,MF_meanfspk,'b','LineWidth',3)
box off
set(gca,'FontSize',16);
ylim([0 maxMF])
xlim([startingzoom endingzoom])
set(h1,'Position',[0.13 0.767258064516129 0.775 0.21688819795310238]);
set(h2,'Position',[0.13 0.548172043010753 0.775 0.21688819795310238]);
set(h3,'Position',[0.13 0.329086021505376 0.775 0.21688819795310238]);
set(h4,'Position',[0.13 0.11 0.775 0.21688819795310238]);
end
end
%% 3D PLOTS
seiv=1;
for starting = 0:5:300
ending = starting+0.5;
output_cut=output(output(:,1)>=starting,:);
output_cut=output_cut(output_cut(:,1)<ending,:);
% PC FREQUENCY ANALYSIS
fspk=zeros((ending-starting)*1000,PCNum);
for i=1:PCNum
for j=1:size(fspk,1)
tspk=output_cut(output_cut(:,2)==i-1+FirstPC,1);
tspk=ceil(tspk.*1000);
inizio=j;
fine=j+window;
a=tspk(tspk>=inizio+starting*1000);
a=a(a<fine+starting*1000);
fspk(j,i)=numel(a)*1000/window;
clear a
end
end
PC_meanfspk=smooth(mean(fspk,2),mobile);
PC_save(seiv,:)=PC_meanfspk;
clear fspk
% DCN FREQUENCY ANALYSIS
fspk=zeros((ending-starting)*1000,DCNNum);
for i=1:DCNNum
for j=1:size(fspk,1)
tspk=output_cut(output_cut(:,2)==i-1+FirstDCN,1);
tspk=ceil(tspk.*1000);
inizio=j;
fine=j+window;
a=tspk(tspk>=inizio+starting*1000);
a=a(a<fine+starting*1000);
fspk(j,i)=numel(a)*1000/window;
clear a
end
end
DCN_meanfspk1=smooth(mean(fspk(:,1:DCNNum/2),2),mobile);
DCN_meanfspk2=smooth(mean(fspk(:,DCNNum/2+1:DCNNum),2),mobile);
DCN_meanfspk=DCN_meanfspk1-DCN_meanfspk2;
DCN_save(seiv,:)=DCN_meanfspk;
clear fspk
seiv=seiv+1;
end
figure
set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.
j=1;
for i=1:seiv-1
MatrixDCN(j,1:500)=DCN_save(j,:);
plot3(10*(j-1)*ones(size(10:500)),10:500,MatrixDCN(j,10:500),'k')
hold on
j=j+1;
end
box off
ylabel({'Time[ms]'},'FontSize',16);
xlabel({'Trial'},'FontSize',16);
zlabel({'Firing Rate [Hz]'},'FontSize',16);
ylim([0 500])
xlim([0 600])
set(gca,'FontSize',16);
view([222 57]);
set(gca,'XDir','reverse');
set(gca,'YDir','reverse');
figure
set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.
j=1;
for i=1:seiv-1
MatrixPC(j,1:500)=PC_save(j,:);
plot3(10*(j-1)*ones(size(1:500)),1:500,MatrixPC(j,1:500),'g')
hold on
j=j+1;
end
box off
ylabel({'Time[ms]'},'FontSize',16);
xlabel({'Trial'},'FontSize',16);
zlabel({'Firing Rate [Hz]'},'FontSize',16);
ylim([0 500])
xlim([0 600])
set(gca,'FontSize',16);
view([222 57]);
set(gca,'XDir','reverse');
set(gca,'YDir','reverse');
%% CR PERCENTAGE
count=1;
nrep=600;
CR_ordered(1,1)=CR(1);
for i=2:size(CR,1)
if CR(i)~=CR(i-1)
CR_ordered(count+1,1)=CR(i);
count=count+1;
end
end
count=1;
for i=1:nrep
CR_ordered2=CR_ordered(find(CR_ordered<=(i/2)));
CRS=CR_ordered2(find(CR_ordered2>(i/2-5)))
CR_event(count)=size(CRS,1)*10;
count=count+1;
end
figure
set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.
plot(1:nrep,CR_event,'k','Marker','pentagram','MarkerFaceColor',[0 0 0]);
% Create xlabel
xlabel({'Trial'},'FontSize',15);
% Create ylabel
ylabel({'% of CR'},'FontSize',15);
% Create Title
title({'Percentage of Conditioned Responses (Computed every trial)'},'LineWidth',3,'FontSize',19,...
'FontName','Agency FB');
set(gca,'FontSize',15);
%% SYNAPTIC WEIGHTS
cellnumber=250;
for i=5:5:300
filename=sprintf('FinalWeights%d.dat',i);
load(filename);
end
FinalWeights0=load('../weights_EBCC.dat');
total = 0;
for i=1:size(FinalWeights10,1)
total=total+FinalWeights10(i,1);
end
pesi=zeros(total,60);
count = 1;
j=0;
eval([ 'M = FinalWeights' num2str(j) ';' ]);
for i=1:size(M,1)
if M(i,1)==1
pesi(count,1)=M(i,2);
count = count+1;
else
while(M(i,1)>=1)
pesi(count,1)=M(i,2);
M(i,1) = M(i,1)-1;
count = count+1;
end
end
end
for j=0.5:0.5:30
count = 1;
eval([ 'M = FinalWeights' num2str(j*10) ';' ]);
for i=1:size(M,1)
if M(i,1)==1
pesi(count,j*2+1)=M(i,2);
count = count+1;
else
while(M(i,1)>=1)
pesi(count,j*2+1)=M(i,2);
M(i,1) = M(i,1)-1;
count = count+1;
end
end
end
end
grey = [0.45,0.45,0.45];
figure
set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.
for i=5263:10:34132
if(pesi(i,nrep/50)>=pesi(i,1))
plot([0:10:nrep],pesi(i,:),'r')
else
plot([0:10:nrep],pesi(i,:),'Color',grey)
end
s = ['Progress: ' num2str(i/total*100) '%'];
hold on
display(s)
end
xlabel({'Trial'},'FontSize',15);
% Create ylabel
ylabel({'GR-PC Synaptic Weights [nS]'},'FontSize',15);
% Create Title
title({'GR-PC Synaptic Weights'' Evolution'},'LineWidth',3,'FontSize',19,...
'FontName','Agency FB');
set(gca,'FontSize',15)
xlim([0 nrep])