% 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 Supplement Fig. S7 (drift-diffusion matched 
%synaptic distributions)

clear variables;
th=20*10^-3;%V Threshold
% lk=1/(20*10^-3);%s^-1 leak constant 1/RC
% dt=0.2*10^-3;%s time step
% tf=200*10^-3; %s final time
dv0=0.2*10^-3;%V defines max step in potential
v0=(0:dv0:th);%V defines the grid for the membrane potential
% p0=zeros(size(v0));
% p0(1)=1;      
% tvec=0:dt:tf;
% voltage_fraction=1;


fnt = 18;
fntlab = 10;
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');


%%%%%ROW 1: matched distributions with panels for synaptic weight plot,
%%%%%equilibrium distributions, and firing rates (matching of both drift and
%%%%%diffusion)
%%%To be placed in the legend:
%%%Firing rates:
%%%Gauss =      745.9868 Hz
%%%Exp =        952.3808 Hz
%%%Hyperbolic = 4166.7205 Hz
%%%Lognormal =  1062.1518 Hz
%%%Delta =      526.3158 Hz
%%%Bimodal =    1825.33 Hz
%%%All means =  0.001 V (all drifts = 1 V/s)
%%%All diffusionss: (=firing rate * sum(probabilites*weights^2)) = 0.0019
%%%V^2/s

match_drift_only=0;
match_4_distributions_1_var;
qdw_gauss=[wts_gauss;ones(size(wts_gauss));qwts_gauss];
qdw_exp=[wts_exp; ones(size(wts_exp)); qwts_exp]; 
qdw_hyp=[wts_hyp; ones(size(wts_hyp)); qwts_hyp]; 
qdw_logn=[wts_logn; ones(size(wts_logn)); qwts_logn]; 
qdw_del=[0.0018; 1; 1];
qdw_bidel=[0.0005 0.015;1 1;0.9967 0.0033];
fr_del=1/qdw_del(1,1);
fr_bidel=1825.33;

%%%%%SYNAPTIC DISTRIBUTION PLOTS
figure;
s1=subplot1(1,3,1);
minyaxis=10^-5;
dplot=wts_hyp(1)/10;
wts_del=qdw_del(1,1)-dplot:dplot:qdw_del(1,1)+dplot;
qwts_del=0.1*minyaxis*ones(size(wts_del));
qwts_del(2)=1;
wts_bidel=[qdw_bidel(1,1)-dplot:dplot:qdw_bidel(1,1)+dplot qdw_bidel(1,2)-dplot:dplot:qdw_bidel(1,2)+dplot];
qwts_bidel=0.1*minyaxis*ones(size(wts_bidel));
qwts_bidel(2)=qdw_bidel(3,1);
qwts_bidel(5)=qdw_bidel(3,2);

semilogy(wts_del,qwts_del,'.-',wts_gauss,qwts_gauss,'.-',wts_exp,qwts_exp,'.-',wts_logn,qwts_logn,'.-','MarkerSize',15);
hold on;
plot(wts_bidel,qwts_bidel,'.-','Color',[0.2 0 0.5],'MarkerSize',15);
plot(wts_hyp,qwts_hyp,'.-','Color',[0.8 0.8 0],'MarkerSize',15);
axis([0 0.02 minyaxis 1.2]);
set(gca,'FontSize',fnt);
set(gca,'XTick',ytc);
set(gca,'XTickLabel',ytl);
set(gca,'FontSize',fnt);
ylabel('Probability');xlabel('synaptic wt (mV) \newline            (a)');
axis square;
leg1 = legend('Delta','Gauss','Exp','LogN','Bimod','PL','Location','North');
set(leg1,'FontSize',fntlab,'Location','North','Interpreter','tex');

load('SuppFigS7_distributions.mat');
load('SuppFigS7_frates.mat');

%%%%%EQUILIBRIUM DISTRIBUTION PLOTS
s2=subplot1(1,3,2);
scale = 10^-3;
plot(v0,scale*prob_del_unif(end,:),v0,scale*prob_gauss_unif(end,:),v0,scale*prob_exp_unif(end,:),v0,scale*prob_logn_unif(end,:),'linewidth',2);
hold on;
plot(v0,scale*prob_bidel_unif(end,:),'color',[0.2,0,0.5],'linewidth',2);
plot(v0,scale*prob_hyp_unif(end,:),'color',[0.8,0.8,0],'linewidth',2);
axis([0 v0(end) 0 scale*200]);
set(gca,'FontSize',fnt);
set(gca,'FontSize',fnt);
ylabel('Prob density (mV^{-1})');xlabel('voltage (mV) \newline         (b)');
axis square;
leg2 = legend('Delta','Gauss','Exp','LogN','Bimod','PL','Location','NorthWest');
set(leg2,'FontSize',fntlab,'Location','NorthWest','Interpreter','tex');


%%%%%FIRING RATE PLOTS
s3=subplot1(1,3,3);
plot(t_vec_del,fout_del,t_vec_gauss,fout_gauss,t_vec_exp,fout_exp,t_vec_logn,fout_logn,'linewidth',2);
hold on;
plot(t_vec_bidel,fout_bidel,'color',[0.2,0,0.5],'linewidth',2);
plot(t_vec_hyp,fout_hyp,'color',[0.8,0.8,0],'linewidth',2);
set(gca,'FontSize',fnt);
set(gca,'XTick',xtc);
set(gca,'XTickLabel',xtl);
set(gca,'FontSize',fnt);
ylabel('Firing rate (Hz)');xlabel('time (ms) \newline      (c)');
axis square;
leg3 = legend('Delta','Gauss','Exp','LogN','Bimod','PL','Location','SouthEast');
set(leg3,'FontSize',fntlab,'Location','SouthEast','Interpreter','tex');

print2eps Fig_S7.eps
