Adaptive robotic control driven by a versatile spiking cerebellar network (Casellato et al. 2014)

 Download zip file 
Help downloading and running models
Accession:167414
" ... We have coupled a realistic cerebellar spiking neural network (SNN) with a real robot and challenged it in multiple diverse sensorimotor tasks. ..."
Reference:
1 . Casellato C, Antonietti A, Garrido JA, Carrillo RR, Luque NR, Ros E, Pedrocchi A, D'Angelo E (2014) Adaptive robotic control driven by a versatile spiking cerebellar network. PLoS One 9:e112265 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Cerebellum;
Cell Type(s): Cerebellum Purkinje cell; Cerebellum interneuron granule cell; Cerebellum golgi cell;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: EDLUT;
Model Concept(s): Detailed Neuronal Models; Synaptic Plasticity; Long-term Synaptic Plasticity; Action Potentials; Learning; STDP;
Implementer(s): D'Angelo, Egidio [dangelo at unipv.it]; Garrido, Jesus A [jesus.garrido at unipv.it]; Luque, Niceto R. [nluque at ugr.es];
Search NeuronDB for information about:  Cerebellum Purkinje cell; Cerebellum interneuron granule cell;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %          RESULT ANALYSIS FOR EBCC SIMULATION ISI = 200 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])