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

clear variables;

%Parameters
th=20*10^-3;%V Threshold
lk=1/(20*10^-3);%s^-1 leak constant 1/RC

dt0=0.2*10^-3;%0.2*10^-3;
tf=200*10^-3; %s final time
tvec=0:dt0:tf;

dv0=0.2*10^-3;%V defines max step in potential
v0=(0:dv0:th);%V defines the grid for the membrane potential

scale = 10^3; % to convert from V to mV

flag_poisson=1;             %Poisson input   
flag_cond_syn=0;            % implies current-based synapses
neur_type = 'LIF';

voltage_fraction=1;

%Initial probability density

%delta function probability distribution starting at 0 mV
p0=zeros(size(v0));
p0(v0==0)=1;

[cellpop connect]=init_cellpop_no_dep;

cellpop(1).v0=v0;
cellpop(2).v0=v0;
cellpop(1).p0=p0;
cellpop(2).p0=p0;
cellpop(2).external_fr=zeros(size(tvec));

connect(1).flag_cond_syn=flag_cond_syn;
connect(1).flag_poisson=flag_poisson;
connect(1).voltage_frac=voltage_fraction;

fnt=10;
figno=1;


%%%%%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.333 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;

cellpop(2).external_fr=zeros(size(tvec));

%1) delta-fn distribution of synaptic weights
cellpop(1).external_fr=[fr_del*ones(size(tvec))];
connect(1).qdw=qdw_del;

[cellpop connect t_vec]=dipde_multi_pop_fig1(cellpop,connect,dt0,tf,neur_type);

t_vec_del = t_vec;
fout_del = cellpop(2).fout;
fout_del_up = cellpop(2).fout_up;
fout_del_low = cellpop(2).fout_low;
prob_del_unif = cellpop(2).prob_unif;

%2) bimodal distribution of synaptic weights
cellpop(1).external_fr=[fr_bidel*ones(size(tvec))];
connect(1).qdw=qdw_bidel;

[cellpop connect t_vec]=dipde_multi_pop_fig1(cellpop,connect,dt0,tf,neur_type);

t_vec_bidel = t_vec;
fout_bidel = cellpop(2).fout;
fout_bidel_up = cellpop(2).fout_up;
fout_bidel_low = cellpop(2).fout_low;
prob_bidel_unif = cellpop(2).prob_unif;

%3) Gaussian distribution of synaptic weights
cellpop(1).external_fr=[fr_gauss*ones(size(tvec))];
connect(1).qdw=qdw_gauss;

[cellpop connect t_vec]=dipde_multi_pop_fig1(cellpop,connect,dt0,tf,neur_type);

t_vec_gauss = t_vec;
fout_gauss = cellpop(2).fout;
fout_gauss_up = cellpop(2).fout_up;
fout_gauss_low = cellpop(2).fout_low;
prob_gauss_unif = cellpop(2).prob_unif;

%4) Exponential distribution of synaptic weights
cellpop(1).external_fr=[fr_exp*ones(size(tvec))];
connect(1).qdw=qdw_exp;

[cellpop connect t_vec]=dipde_multi_pop_fig1(cellpop,connect,dt0,tf,neur_type);

t_vec_exp = t_vec;
fout_exp = cellpop(2).fout;
fout_exp_up = cellpop(2).fout_up;
fout_exp_low = cellpop(2).fout_low;
prob_exp_unif = cellpop(2).prob_unif;

%5) lognormal distribution of synaptic weights
cellpop(1).external_fr=[fr_logn*ones(size(tvec))];
connect(1).qdw=qdw_logn;

[cellpop connect t_vec]=dipde_multi_pop_fig1(cellpop,connect,dt0,tf,neur_type);

t_vec_logn = t_vec;
fout_logn = cellpop(2).fout;
fout_logn_up = cellpop(2).fout_up;
fout_logn_low = cellpop(2).fout_low;
prob_logn_unif = cellpop(2).prob_unif;

%6) hyperbolic distribution of synaptic weights
cellpop(1).external_fr=[fr_hyp*ones(size(tvec))];
connect(1).qdw=qdw_hyp;

[cellpop connect t_vec]=dipde_multi_pop_fig1(cellpop,connect,dt0,tf,neur_type);

t_vec_hyp = t_vec;
fout_hyp = cellpop(2).fout;
fout_hyp_up = cellpop(2).fout_up;
fout_hyp_low = cellpop(2).fout_low;
prob_hyp_unif = cellpop(2).prob_unif;

%Save
save('SuppFigS7_distributions.mat','prob_gauss_unif','prob_exp_unif','prob_hyp_unif','prob_logn_unif','prob_del_unif','prob_bidel_unif');
save('SuppFigS7_frates.mat','t_vec_gauss','fout_gauss','t_vec_exp','fout_exp','t_vec_hyp','fout_hyp','t_vec_logn','fout_logn','t_vec_del','fout_del','t_vec_bidel','fout_bidel');
save('SuppFigS7_frates_up.mat','t_vec_gauss','fout_gauss_up','t_vec_exp','fout_exp_up','t_vec_hyp','fout_hyp_up','t_vec_logn','fout_logn_up','t_vec_del','fout_del_up','t_vec_bidel','fout_bidel_up');
save('SuppFigS7_frates_low.mat','t_vec_gauss','fout_gauss_low','t_vec_exp','fout_exp_low','t_vec_hyp','fout_hyp_low','t_vec_logn','fout_logn_low','t_vec_del','fout_del_low','t_vec_bidel','fout_bidel_low');
