% 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. S6 (matched Gaussian
%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

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;

%%%ROW 1: Synaptic distributions set as two Gaussians with two means and two 
%standard deviations

%Parameter sets to run:
mean_std_matrix=[6 6 3 3;2 1 1 0.5]*10^-3;    %First row: means (V), second row: standard deviations (V)
fr=[250 250 500 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=[];
v0=(0:dv0:th);%V defines the grid for the membrane potential
p0=zeros(size(v0));
[val zero_index]=min(abs(v0));
p0(zero_index)=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;

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(mean_std_matrix(1,ii)) 'V, SD=' num2str(mean_std_matrix(2,ii)) 'V'];
end;

save('SuppFigS6_1.mat','wts_g','qdw_g','fout_g','tvec','fout_g_up','fout_g_low','legendstring');


%%%ROW 2: Synaptic distributions set as two Gaussians with means on either
%%%side of zero and equal standard deviations

%Double Gaussians
%Parameter sets to run:
mean_std_matrix=[6 6 3 3;2 1 1 0.5]*10^-3;    %First row: means (V), second row: standard deviations (V)
num_trials=length(mean_std_matrix(1,:));
fr=500;                     %Hz
wts_g=-th:dv0:th;
qdw_g=zeros(num_trials,length(wts_g));
fout_g=[];
fout_g_up=[];
fout_g_low=[];
equil_g=[];
v0=-5*th:dv0:th;
p0=zeros(size(v0));
[val zero_index]=min(abs(v0));
p0(zero_index)=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;

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*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(mean_std_matrix(1,ii)) 'V, SD=' num2str(mean_std_matrix(2,ii)) 'V'];
end;

save('SuppFigS6_2.mat','wts_g','qdw_g','fout_g','tvec','fout_g_up','fout_g_low','legendstring');
