% 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 Figure 2 plots (matched synaptic distributions)

clear variables;

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=100*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;


%%%ROW 1: Synaptic distributions set as two self-similar Gaussians with two means
%Parameter sets to run:
mean_std_matrix=[6 3;2 1]*10^-3;    %First row: means (V), second row: standard deviations (V)
fr=[250 500];   %Hz
num_trials=length(mean_std_matrix(1,:));
wts_g=0:dv0:th;
qdw_g=zeros(num_trials,length(wts_g));
fout_g=[];
fout_g_up = [];
fout_g_low = [];
equil_g=[];
legendstring=cell(num_trials,1);
 for ii=1:num_trials
    qdw_g(ii,:)=exp(-(mean_std_matrix(1,ii)-wts_g).^2/(2*mean_std_matrix(2,ii)^2))+exp(-(-mean_std_matrix(1,ii)-wts_g).^2/(2*mean_std_matrix(2,ii)^2));
    qdw_g(ii,:)=qdw_g(ii,:)/sum(qdw_g(ii,:));
    qdw_temp=[wts_g;ones(size(wts_g));qdw_g(ii,:)];
    cellpop(1).external_fr=[fr(ii)*ones(size(tvec))];
    connect(1).qdw=qdw_temp;
    [cellpop connect t_vec]=dipde_multi_pop_fig1(cellpop,connect,dt0,tf,neur_type);
    fout_g=[fout_g; cellpop(2).fout];
    fout_g_up = [fout_g_up;cellpop(2).fout_up];
    fout_g_low = [fout_g_low;cellpop(2).fout_low];
    equil_g=[equil_g; cellpop(2).prob_unif(end,:)];
    legendstring{ii}=['Mean=' num2str(scale*mean_std_matrix(1,ii)) 'mV,\newline SD=' num2str(scale*mean_std_matrix(2,ii)) 'mV,\newline f_{in} =' num2str(fr(ii)) 'Hz'];
 end;

 
%%%ROW 1: Synaptic distributions set as two self-similar Gaussians with two
%%%means with LOWER frequencies
%Parameter sets to run:
mean_std_matrix2=[6 3;2 1]*10^-3;    %First row: means (V), second row: standard deviations (V)
fr2=[100 200];   %Hz
num_trials2=length(mean_std_matrix2(1,:));
wts_g2=0:dv0:th;
qdw_g2=zeros(num_trials2,length(wts_g2));
fout_g2=[];
fout_g2_up = [];
fout_g2_low = [];
equil_g2=[];
legendstring2=cell(num_trials2,1);
 for ii=1:num_trials2
    qdw_g2(ii,:)=exp(-(mean_std_matrix2(1,ii)-wts_g2).^2/(2*mean_std_matrix2(2,ii)^2))+exp(-(-mean_std_matrix2(1,ii)-wts_g2).^2/(2*mean_std_matrix2(2,ii)^2));
    qdw_g2(ii,:)=qdw_g2(ii,:)/sum(qdw_g2(ii,:));
    qdw_temp2=[wts_g2;ones(size(wts_g2));qdw_g2(ii,:)];
    cellpop(1).external_fr=[fr2(ii)*ones(size(tvec))];
    connect(1).qdw=qdw_temp2;
    [cellpop connect t_vec]=dipde_multi_pop_fig1(cellpop,connect,dt0,tf,neur_type);
    fout_g2=[fout_g2; cellpop(2).fout];
    fout_g2_up = [fout_g2_up; cellpop(2).fout_up];
    fout_g2_low = [fout_g2_low; cellpop(2).fout_low];
    equil_g2=[equil_g2; cellpop(2).prob_unif(end,:)];
    legendstring2{ii}=['Mean=' num2str(scale*mean_std_matrix2(1,ii)) 'mV,\newline SD=' num2str(scale*mean_std_matrix2(2,ii)) 'mV,\newline f_{in} =' num2str(fr(ii)) 'Hz'];
end;

%--------------------------------------------------------------------

%%ROW 2: matched distributions with panels for synaptic weight plot,
%%equilibrium distributions, and firing rates

fr = 1000;

match_drift_only=1;
match_4_distributions_1_var;
qdw_del=[0.001; 1; 1];
qdw_bidel=[0.0005 0.015;1 1;14/14.5 1-(14/14.5)];
qdw_hyp=[wts_hyp; ones(size(wts_hyp)); qwts_hyp]; 
qdw_gauss=[wts_gauss;ones(size(wts_gauss));qwts_gauss];
qdw_exp=[wts_exp; ones(size(wts_exp)); qwts_exp]; 
qdw_logn=[wts_logn; ones(size(wts_logn)); qwts_logn]; 

cellpop(1).external_fr=[fr*ones(size(tvec))];
cellpop(2).external_fr=zeros(size(tvec));

%1) delta-fn distribution of synaptic weights
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
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
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
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
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
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('Fig2_gaussian_plots.mat','v0','mean_std_matrix','wts_g','qdw_g','equil_g','equil_g2','tvec','fout_g','fout_g_up','fout_g_low','fout_g2','fout_g2_up','fout_g2_low');
save('Fig2_matched_distributions.mat','prob_gauss_unif','prob_exp_unif','prob_hyp_unif','prob_logn_unif','prob_del_unif','prob_bidel_unif');
save('Fig2_matched_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('Fig2_matched_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('Fig2_matched_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');

