% 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. S14 (effects of external fluctuations for a given
%%%%set of distributions with input rate of 8000 Hz). Please refer to 
%METHODS: Fluctuation Analysis in the paper for the implementation.

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

dt=0.5*10^-3;%s time step
tf=300*10^-3; %s final time
tvec=0:dt:tf;

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;      

voltage_fraction=1;

num_inputs=100;     %number of fluctuations
proportion=0.1;     %proportion of neurons receiving input for a 'single' fluctuation

%%%Get distributions
match_drift_only=1;
match_4_distributions_1_var;
qdw_gauss=[wts_gauss;ones(size(wts_gauss));qwts_gauss];
qdw_hyp=[wts_hyp; ones(size(wts_hyp)); qwts_hyp]; 
qdw_exp=[wts_exp; ones(size(wts_exp)); qwts_exp]; 
qdw_logn=[wts_logn; ones(size(wts_logn)); qwts_logn]; 
qdw_del=[0.001; 1; 1];
qdw_bidel=[0.0005 0.015;1 1;14/14.5 1-(14/14.5)];
if match_drift_only~=1
    qdw_del=[0.0018; 1; 1];
    qdw_bidel=[0.0005 0.015;1 1;0.9967 0.0033];
end
qdw(1).mat=qdw_gauss;
qdw(2).mat=qdw_hyp;
qdw(3).mat=qdw_exp;
qdw(4).mat=qdw_logn;
qdw(5).mat=qdw_del;
qdw(6).mat=qdw_bidel;
fr=fr_gauss*ones(length(qdw),1);
if match_drift_only~=1
    fr=[fr_gauss fr_hyp fr_exp fr_logn 1/qdw_del(1,1) 1825.3];
end;
fr=fr*8;   %input rate becomes 8000 Hz

%%%%%ROW 3: Generate relative excitability plots for all distributions with
%%%%%and without depression
excit_no_sd=zeros(num_inputs,6);
excit_sd=zeros(num_inputs,6);
savefiles=cell(6,1);
savefiles{1}='gauss_fluct_8000.mat';
savefiles{2}='hyper_fluct_8000.mat';
savefiles{3}='exp_fluct_8000.mat';
savefiles{4}='logn_fluct_8000.mat';
savefiles{5}='del_fluct_8000.mat';
savefiles{6}='bidel_fluct_8000.mat';

for dist_type=1:6
    
    %Run with no depression until equilibrium
    [cellpop connect]=init_cellpop_no_dep;
    
    connect(1).qdw=qdw(dist_type).mat;
    cellpop(1).external_fr=fr(dist_type)*ones(size(tvec));
    
    [cellpop connect tvec]=dipde_multi_pop2(cellpop,connect,dt,tf,1);

    %fluctuations with varying numbers of inputs (no leak, no depression)
    dt=tvec(2)-tvec(1);
    
    [newcellpop newconnect]=init_cellpop_no_dep;
    
    newcellpop(1).external_fr=(proportion/dt)*ones(num_inputs+1,1);
    newconnect(1).flag_poisson=0;
    newconnect(1).qdw=qdw(dist_type).mat;
    newcellpop(2).p0=cellpop(2).prob_unif(end,:);
    
    [newcellpop newconnect tvec2]=dipde_multi_pop2(newcellpop,newconnect,dt,dt*num_inputs,0);   
    excit_no_sd(:,dist_type)=newcellpop(2).fout(2:end)'/newcellpop(2).fout(2);

    %fluctuations with varying numbers of inputs (no leak, with depression)
    dt=tvec(2)-tvec(1);
    
    [newcellpop_sd newconnect_sd]=init_cellpop_no_dep;
    
    newcellpop_sd(1).external_fr=(proportion/dt)*ones(num_inputs+1,1);
    newcellpop_sd(1).u_pts=0:0.01:1;
    [maxval maxind]=max(newcellpop(2).fout);
    newcellpop_sd(1).u_frac=((newcellpop(2).fout(2)/maxval)^(1/(maxind*proportion)));
    newcellpop_sd(1).u_tau=100*10^-3;
    newconnect_sd(1).flag_poisson=0;
    newconnect_sd(1).nsyn=1;
    newconnect_sd(1).qdw=qdw(dist_type).mat;
    newcellpop_sd(2).p0=cellpop(2).prob_unif(end,:);
    
    [newcellpop_sd newconnect_sd tvec2]=dipde_multi_pop2(newcellpop_sd,newconnect_sd,dt,dt*num_inputs,0);   
    excit_sd(:,dist_type)=newcellpop_sd(2).fout(2:end)'/newcellpop_sd(2).fout(2);
    
    save(savefiles{dist_type},'newcellpop','newcellpop_sd');
end;

%save('excitabilities_8000.mat','excit_no_sd','excit_sd');
