% This file is licensed under the Allen Institute Terms of Use, which can be found at:
%http://www.alleninstitute.org/Media/policies/terms_of_use_content.html
%------------------------------------------------------------------------------------

%Script to generate Supplement Fig. S1 showing convergence of 3 methods for 
%1mV and 1000Hz input for current-based synapses

clear variables;

% Settings for plot
fnt = 18;
xtc = [0;50;100;150;200]*10^-3;  %x-axis ticks for time in (ms)
ytc = [0;5;10;15;20]*10^-3; %y-axis ticks for membrane potential in (mV)
xtl = [0;50;100;150;200]; %x-axis labels for time in (ms)
ytl = [0;5;10;15;20]; %y-axis labels for membrane potential in (mV)

set(0,'DefaultFigureWindowStyle','docked');
set(0,'defaultAxesFontName', 'Arial');

load 'SuppFigS1.mat';

figure;

%%%Probability distribution time-evolution plots

%Time-evolution of probability distribution of membrane potential from NEST simulations
s1 = subplot1(2,3,1);
imagesc(tsim1,vsim1,probsim_curr);set(gca,'YDir','normal');shading flat;
set(gca, 'FontSize',fnt); 
set(gca,'XTick',xtc); set(gca,'YTick',ytc);
set(gca,'XTickLabel',xtl); set(gca,'YTickLabel',ytl);
set(gca,'FontSize',fnt);
caxis([0 100]);
ylabel('voltage (mV)');
xlabel('time (ms) \newline    (a)');
axis square tight;

%Time evolution of probability distribution of membrane potential from DiPDE
s2 = subplot1(2,3,2);
imagesc(t_curr_dipde,vsim1,prob_unif_curr_dipde(:,1:end-1)');set(gca,'YDir','normal');shading flat;
set(gca,'FontSize',fnt); 
set(gca,'XTick',xtc); set(gca,'YTick',ytc);
set(gca,'XTickLabel',xtl); set(gca,'YTickLabel',ytl);
set(gca,'FontSize',fnt);
caxis([0 100]);
xlabel('time (ms) \newline    (b)');
axis square tight;

%Time evolution of probability distribution of membrane potential from Fokker-Planck
s3 = subplot1(2,3,3);
imagesc(t_fp,vsim1,prob_unif_fp0(:,1:end-1)');set(gca,'YDir','normal');shading flat;
set(gca,'FontSize',fnt); 
set(gca,'XTick',xtc); set(gca,'YTick',ytc);
set(gca,'XTickLabel',xtl); set(gca,'YTickLabel',ytl);
set(gca,'FontSize',fnt);
caxis([0 100]);
xlabel('time (ms) \newline    (c)');
colorbar('location','EastOutside','YTickLabel',{'0','0.2','0.4','0.6','0.8','1'});
set(gca,'FontSize',fnt);
axis square tight;

%Settings to align heat-maps properly
s1Pos = get(s1,'position');
s2Pos = get(s2,'position');
s3Pos = get(s3,'position');
s2Pos(2) = s1Pos(2) + 0.5*0.0787;
s3Pos(2) = s1Pos(2) + 0.5*0.0787;
s2Pos(3:4) = [s1Pos(3:4)];
s3Pos(3:4) = [s1Pos(3:4)];
set(s2,'position',s2Pos);
set(s3,'position',s3Pos);


%FIRING RATE PLOTs
%-----------------
s4 = subplot1(2,3,4);
plot(tplot, simfrate_curr(1:end-1),'b','LineWidth',2);
set(gca,'XTick',xtc);
set(gca,'XTickLabel',xtl);
hold on;
plot(t_curr_dipde,fout_curr_dipde, '--r', 'LineWidth',2);
plot(t_fp,fout_fp/dt_fp, '--', 'LineWidth',2,'Color',[0,0.5,0]);
legend('Sim','DiPDE','FP','Location','SouthEast');
set(gca, 'FontSize',fnt);
set(gca,'FontSize',fnt);
ylabel('Firing rate (Hz)');xlabel('time (ms) \newline    (d)');
axis tight;

%Membrane potential probability distribution at t = 200 ms
%---------------------------------------------------------
probsim_curr(end+1,end-1) = 0;
s5 = subplot1(2,3,5);
plot(v0(1:end),10^-3*probsim_curr(1:end,end-1),'b', 'LineWidth',2);
hold on;
plot(v0(1:end),10^-3*prob_unif_curr_dipde(end-1,1:end),'--r', 'LineWidth',2);
plot(v0(1:end),10^-3*prob_unif_fp0(end-1,1:end),'--', 'LineWidth',2,'Color',[0,0.5,0]);
set(gca,'XTick',ytc);
set(gca,'XTickLabel',ytl);
legend('Sim','DiPDE','FP','Location','NorthWest');
set(gca, 'FontSize',fnt);
ylabel('Prob density (mV^{-1})');xlabel('voltage (mV) \newline        (e)');
axis tight;

%Align panels S1(d) and S1(e) properly
%-----------------------------------
s4Pos = get(s4,'position');
s5Pos = get(s5,'position');
s4Pos(1) = s1Pos(1);
s5Pos(1) = s2Pos(1) + 0.8*s2Pos(3);
s4Pos(3:4) = [1.5*s1Pos(3) s1Pos(4)];
s5Pos(3:4) = [1.5*s2Pos(3) s2Pos(4)];
s4Pos(3:4) = s5Pos(3:4);
set(s4,'position',s4Pos);
set(s5,'position',s5Pos);


print2eps Fig_S1.eps
