% 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 plots for Figure 3 (input-output curves for a given
%%%%set of distributions.

clear variables;

baserates=0.5:0.5:2;    %Set of input rates to plot (will be multiplied by firing rate)
fractions=0.5:0.1:1;    %Set of fractions of firing rates at which to take equilibrium distributions to identify peak transients
baserate=1000;

%Plot settings
fnt = 18;
fntlab = 10;

xtc = [0;500;1000]*10^-3;
xtl = [0;500;1000];

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


%%%%%PANEL A: Description of simulations

figure;
s1=subplot1(2,3,1);
legendstring=cell(6,1);
hold all;
for ii=1:length(fractions)-2
    plot([0 500 500 1000],[fractions(ii)*baserate fractions(ii)*baserate baserate baserate],'linewidth',2);
    legendstring{ii}=['Frac = ' num2str(fractions(ii))];
end;
plot([0 500 500 1000],[fractions(end-1)*baserate fractions(end-1)*baserate baserate baserate],'Color',[0.2,0,0.5],'linewidth',2);
    legendstring{5}=['Frac = ' num2str(fractions(5))];
plot([0 1000], [baserate,baserate],'Color',[0.8,0.8,0],'linewidth',2);
axis([0 1000 0 1500])
hold off;
legendstring{length(fractions)}=['Frac = ' num2str(fractions(end))];
set(gca,'FontSize',fnt);
set(gca,'FontSize',fnt);
ylabel('Input firing rate (Hz)');xlabel('time (ms) \newline    (a)');
axis square; box on;
leg1 = legend(legendstring,'Location','Southeast');
set(leg1,'FontSize',fntlab,'Interpreter','tex');

%%%PANEL B: Sample output firing rates for delta function distribution
load('delta_all_firing.mat');
s2 = subplot1(2,3,2);
hold all;
for ii=1:length(fractions)-2
    plot(tvec_store{ii,3},firing_rate_store{ii,3},'linewidth',2);
    legendstring{ii}=['Frac = ' num2str(fractions(ii))];
end;
plot(tvec_store{5,3},firing_rate_store{5,3},'Color',[0.2,0,0.5],'linewidth',2);
legendstring{5}=['Frac = ' num2str(fractions(5))];
plot(tvec_store{6,3},firing_rate_store{6,3},'Color',[0.8,0.8,0],'linewidth',2);
legendstring{6}=['Frac = ' num2str(fractions(6))];
hold off;
set(gca,'FontSize',fnt);
set(gca,'XTick',xtc);
set(gca,'XTickLabel',xtl);
set(gca,'FontSize',fnt);
ylabel('Output firing rate (Hz)');xlabel('time (ms) \newline    (b)');
axis square; box on;
leg2 = legend(legendstring,'Location','Southeast');
set(leg2,'FontSize',fntlab,'Interpreter','tex');

%%%%%%%%%%PANEL C: Zooming in around jump in input firing rate
load('delta_all_firing.mat');
s3 = subplot1(2,3,3);
hold all;
for ii=1:length(fractions)-2
    plot(tvec_store{ii,3}(490:570),firing_rate_store{ii,3}(490:570),'linewidth',2);
    legendstring{ii}=['Frac = ' num2str(fractions(ii))];
end;
plot(tvec_store{5,3}(490:570),firing_rate_store{5,3}(490:570),'Color',[0.2,0,0.5],'linewidth',2);
legendstring{5}=['Frac = ' num2str(fractions(5))];
plot(tvec_store{6,3}(490:570),firing_rate_store{6,3}(490:570),'Color',[0.8,0.8,0],'linewidth',2);
legendstring{6}=['Frac = ' num2str(fractions(6))];
hold off;
set(gca,'FontSize',fnt);
axis([tvec_store{1,3}(490) tvec_store{1,3}(570) 0 25]);
xtcz = [500;520;540;560]*10^-3;
xtlz = [500;520;540;560];
set(gca,'XTick',xtcz);
set(gca,'XTickLabel',xtlz);
set(gca,'FontSize',fnt);
ylabel('Output firing rate (Hz)');xlabel('time (ms) \newline    (c)');
axis square; 
box on;
leg3 = legend(legendstring,'Location','Southeast');
set(leg3,'FontSize',fntlab,'Interpreter','tex');

%%%%%PANEL D: Peak output rates for different fractions of base input rate with delta function 
load('delta_transient_frates.mat');
s4=subplot1(2,3,4);
hold all;
legendstring=cell(length(fractions),1);
for ii=1:length(fractions)-2
    indices=find(transient_frates(:,1)==fractions(ii));
    plot(transient_frates(indices,2),transient_frates(indices,3),'-o','linewidth',2);
    legendstring{ii}=['Frac = ' num2str(fractions(ii))];
end;
indices=find(transient_frates(:,1)==fractions(5));
plot(transient_frates(indices,2),transient_frates(indices,3),'-o','Color',[0.2,0,0.5],'linewidth',2);
legendstring{5}=['Frac = ' num2str(fractions(5))];
indices=find(transient_frates(:,1)==fractions(6));
plot(transient_frates(indices,2),transient_frates(indices,3),'-o','Color',[0.8,0.8,0],'linewidth',2);
legendstring{6}=['Frac = ' num2str(fractions(6))];
hold off;
set(gca,'FontSize',fnt);
set(gca,'FontSize',fnt);
ylabel('Peak output FR (Hz)');xlabel('Syn input events per sec \newline                 (d)');
axis square tight;box on;
leg4 = legend(legendstring,'Location','NorthWest');
set(leg4,'FontSize',fntlab,'Interpreter','tex');

%%%%PANEL E: Summary of input-output curves for all distributions
filenames=cell(6,1);
filenames{1}='delta_transient_frates.mat';
filenames{2}='gauss_transient_frates.mat';
filenames{3}='exp_transient_frates.mat';
filenames{4}='logn_transient_frates.mat';
filenames{5}='bidel_transient_frates.mat';
filenames{6}='hyp_transient_frates.mat';

s5 = subplot1(2,3,5);
hold all;
for ii=1:length(filenames)-2
    load(filenames{ii});
    %find slopes at each point from fractions
    slopevals=zeros(length(baserates),1);
    for jj=1:length(baserates)
        y1=transient_frates(transient_frates(:,1)==fractions(1) & transient_frates(:,2)==baserates(jj)*baserate,3);
        y2=transient_frates(transient_frates(:,1)==fractions(end) & transient_frates(:,2)==baserates(jj)*baserate,3);
        x1=baserates(jj)*baserate*fractions(1);
        x2=baserates(jj)*baserate*fractions(2);
        slopevals(jj)=(y1-y2)/(x2-x1);
    end;
    plot(baserate*baserates,slopevals,'-o','linewidth',2);
end;
for ii=5
    load(filenames{ii});
    %find slopes at each point from fractions
    slopevals=zeros(length(baserates),1);
    for jj=1:length(baserates)
        y1=transient_frates(transient_frates(:,1)==fractions(1) & transient_frates(:,2)==baserates(jj)*baserate,3);
        y2=transient_frates(transient_frates(:,1)==fractions(end) & transient_frates(:,2)==baserates(jj)*baserate,3);
        x1=baserates(jj)*baserate*fractions(1);
        x2=baserates(jj)*baserate*fractions(2);
        slopevals(jj)=(y1-y2)/(x2-x1);
    end;
    plot(baserate*baserates,slopevals,'-o','Color',[0.2,0,0.5],'linewidth',2);
end;
for ii=6
    load(filenames{ii});
    %find slopes at each point from fractions
    slopevals=zeros(length(baserates),1);
    for jj=1:length(baserates)
        y1=transient_frates(transient_frates(:,1)==fractions(1) & transient_frates(:,2)==baserates(jj)*baserate,3);
        y2=transient_frates(transient_frates(:,1)==fractions(end) & transient_frates(:,2)==baserates(jj)*baserate,3);
        x1=baserates(jj)*baserate*fractions(1);
        x2=baserates(jj)*baserate*fractions(2);
        slopevals(jj)=(y1-y2)/(x2-x1);
    end;
    plot(baserate*baserates,slopevals,'-o','Color',[0.8,0.8,0],'linewidth',2);
end;
hold off;
set(gca,'FontSize',fnt);
set(gca,'FontSize',fnt);
xlim([500 2000]);
ylim([-0.05 0.15]);
ylabel('Norm. diff. in peak output FR')
xlabel('Syn input events per sec \newline                 (e)');
axis square; box on;
leg5 = legend('Delta','Gauss','Exp','LogN','Bimod','PL','Location','Northwest');
set(leg5,'FontSize',fntlab,'Interpreter','tex');

%%%%PANEL F: Input-output curves with synaptic depression
load('depression_rates_025_1msdt.mat');
minxaxis = 10^2;
maxyaxis = 50;
xtclog = [10^2;10^3;10^4;10^5];
xtllog = num2str(xtclog);
s6=subplot1(2,3,6);
set(gca,'xscale','log');
hold all;
legendstring={'Delta','Gauss','Exp','LogN','Bimod','PL'};
for ii=[0 2 4 5]
    indices=find(transient_frates_dep(:,1)==ii);
    semilogx(transient_frates_dep(indices(4:end),2),transient_frates_dep(indices(4:end),3),'-o','linewidth',2);
end
for ii=[1]
    indices=find(transient_frates_dep(:,1)==ii);
    semilogx(transient_frates_dep(indices(4:end),2),transient_frates_dep(indices(4:end),3),'-o','Color',[0.2,0,0.5],'linewidth',2);
end
for ii=[3]
    indices=find(transient_frates_dep(:,1)==ii);
    semilogx(transient_frates_dep(indices(4:end),2),transient_frates_dep(indices(4:end),3),'-o','Color',[0.8,0.8,0],'linewidth',2);
end
hold off;
set(gca,'FontSize',fnt);
set(gca,'XTick',xtclog);
set(gca,'XTickLabel',xtllog);
set(gca,'FontSize',fnt);
ylabel('Equil. output FR (Hz)');
xlabel('Syn input events per sec \newline                  (f)');
axis square tight;box on;
leg6 = legend(legendstring,'Location','NorthWest');
set(leg6,'FontSize',fntlab,'Interpreter','tex');

print2eps Fig3.eps
